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ABSTRACT 

We describe an extensive suite of numerical calculations for the collisional evo¬ 
lution of irregular satellite swarms around 1-300 M e planets orbiting at 120 AU 
in the Fomalhaut system. For 10-100 M 0 planets, swarms with initial masses 
of roughly 1% of the planet mass have cross-sectional areas comparable to the 
observed cross-sectional area of Fomalhaut b. Among 30-300 M ffi planets, our 
calculations yield optically thick swarms of satellites for ages of 1-10 Myr. Obser¬ 
vations with HST and ground-based AO instruments can constrain the frequency 
of these systems around stars in the f3 Pic moving group and possibly other nearby 
associations of young stars. 

Subject headings: planetary systems - planets and satellites: formation - proto¬ 
planetary disks - stars: formation - zodiacal dust - circumstcllar matter 


1. INTRODUCTION 


Fomalhaut b is a planet cand idate on an eccentric orbit at a distance of ~ 120 AU from 


the A-type star Fomalhaut (e.g.. iKalas et al.l 120081: ICurrie et al.l 12012c iGalicher et al.l 120131: 


suggest emission from a cloud of dust instead of a planetary photosphere (e.g.. 

Marengo et al. 

2009; 

Janson et al. 

2012; 

Currie et al. 

2012; 

Kalas et al. 

2013; 

Janson et al. 

201 hi). The ob- 


served level of optical emission r equires grains with a total cross-sectional area of roughl y 


10 23 cm 2 (e.g., 


Kalas et al.l 120081: ICurrie et al.l 12012c IGalicher et al.l l2013t IKalas et al.l 120131) . 
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Two types of collision models can produce a clump of dust emission at large dis¬ 
tances from an A-type star. In the simplest picture, two objects with radii of roughly 


Kalas et al. 

2008; 

Galicher et al. 

2013; 

Kenvon et al. 

2014; 

Lawler et al. 

2015) 


panding at t he escape velocity of a pair of 100 km objects are unresolved on 50-100 yr time 
scales (e.g., iGalicher et al.l l2013t iKenvon et al.l 12014 : lLawler et al.l I2015T ) . Smaller ejected 


particles with larger optic al depth generally expand more rapidly (e.g., iGault et al.l Il963 


Housen fc Holsappld 120031) . If these particles contain a reasonable amount of mass, HST 


or JWST images s hould resolve Fomalhaut b within the next decade (e.g., iTa/mavol 12013 


Kenvon et al.l 120141) . When expanding clouds have an internal velocity dispersion, d ifferen- 
tial motion shears the cloud into a ring (IKenvon fc BromlevH2005t IKenvon et al.ll20141) . Over 
the next decade, HST or JWST observations can also test this prediction. 

Alternative models posit a dis k-shaped or a semi-spherical swarm of irregular satellites 


orbit ing a super-Earth mass planet (IKalas et al.ll2008l: Kennedy fe Wvattll201ll; IKenvon et al. 


2014b In this approach, a collisional cascade among particles with radii r < 500-1000 km 
maintains a large population of dust grains over the 200-400 Myr age of Fomalhaut. Analytic 
results for the long-term evolution favor particles with a power-law size distribution in swarms 
with masses of roughly 0.1 M® around planets with masses M p & 10 M®. However, the 
radius of the largest object in the cascade ( r max ) depends on the slope q of the power-law 
size distribution, where smaller q requires larger r max . Large r max requires very massive 
satellite swarms. 


Observational tests of this model are possible but more complicated (e.g.. Ke n yon et al. 


2 01 4). Radiation pressure from the central star prevents small particles with r < 100 /im 
from remaining bound to the planet. Ejection of these particles produces a distinct trail 
along the planet’s orbit. Simple estima tes suggest the density of small particles ejected from 
massive planets is detectable (IKenvon et al. 2014). Testing this aspect of the model requires 
more detailed analyses of the dust ejected from the satellite swarm. 


To explore this picture in more detail, we examine a suite of numerical simulations for 
spherical swarms of satellites orbiting super-Earth mass planets. Aside from deriving the 
evolution of dust clouds as a function of the mass of the planet and the surrounding satellite 
system, we consider how the amount of dust orbiting the planet depends on the initial radius 
of the largest satellite, the bulk strength of satellites, and the recipe for distributing the debris 
from a collision into lower mass objects. For a standard model, swarms with initial masses 
Md = 0.01 M p and r max « 200-400 km orbiting planets with masses M p = 10-100 M® match 
the observed cross-sectional area of Fomalhaut b. Calculations with weaker satellites allow 
lower mass swarms to match the data. 
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For planets with M p ps 30-300 M®, a ~ 100 AU, and ages of 1-10 Myr, satellite swarms 
have large optical depth r ~ 0.1-1. Relative to the central star, predicted contrast ratios of 


Biller et al. 

2013 

) or HST (e.g., 

Schneider et al. 

2014) 


with ground-based AO systems (e.g. 
can place limits on the frequency of optically thick satellite swarms around planets orbiting 
nearby young stars. 


To connect our results with previous analytic work, we begin our discussion with the 
derivation of a simple analytic model (§5]). After summarizing the numerical approach (§2]), 
we describe the outcomes of simulations as a function of various input parameters O- The 
paper concludes with a brief discussion (j|5]) and a summary of the major results (fj6]) . 


ANALYTIC MODEL 


To interpret observations of debris disks c. 2000. IWvatt fe Dentl (1200211 a/nd lDominik fc Decin 


(2003) developed an analytic m odel for_the long-te rm evolution of a swarm of large solid ob¬ 
jects in a circumstellar disk. Kennedy fe Wvattl (120111) later extended this appro ach to 
spherical swarms of satellites orbiting a massive planet (see also Kennedy et al. 2011). In a 
swarm of satellites, objects with radius r, mass m, and mass density p orbit within a spher¬ 
ical shell with width 5a centered at a distance a from a planet with mass M p and radius R p . 
The planet orbits with semimajor axis a p from a central star with mass M* and luminosity 
L*. Destructive collisions between satellites produce a collisional cascade which slowly grinds 
solids into smaller and smaller objects. Defining an upper mass limit m max for solids partic¬ 
ipating in the cascade, the analytic model yields a simple formula for N max (t), the number 
of these large objects as a function of time. If radiation pressure sets m min , a lower mass 
limit for solids with stable orbits around the pl anet, then the ca scade produ ces a powe r-law 
size distribution between m™ n an d in max (e.g., 


Pohnanvilll969l: IWilliams fe Wetherilllll994l: 


O’Brien fe Greenberdl2003l: Kobavashi fe Tanaka] 120101) . Setting the slope of this size distri¬ 
bution yields another simple formula for the time evolution of the surface area of the dust 
cloud Adit). Adopting optical properties for solids in the cloud yields the dust luminosity 
L d (t). 


2.1. Time evolution 


To deri ve expressions for Nm nx (t ) and A d (t), IWvatt fe Dentl (120021) . iDominik &; Decin 


(120031 ). and Kennedy & Wya t t (20111) adopt the particle-in-a-box model, where kinetic the- 
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ory sets the collision rate N. Defining V as the volume of the spherical shell, a as the 
geometric cross-section, and v as the relative particle velocity, each particle has a collision 
time t c ~ Nq 1 (V/2av) where N 0 is the initial number of particle^]. If all collisions between 
particles are destructive, the number of large particles declines at a rate N ma x ~ ~^ma X /-^otc- 
Solving for N max (t): 


A max (t) 


N, 


max, 0 


where N„ 


1 + t/t c 

C'O is the number of large objects at t — 0. 


( 1 ) 


With N max (t) known, the total cross-sectional area and dust luminosity follow. For any 
power-law size distribution with N{r) oc r~ q , the total mass Md and cross-sectional area Ad 
of the swarm are simple functions of the minimum size, the maximum size, and the slope 
q (e.g., IWvatt fe Dentll2002l; iDominik fe Decinl 1200311 . The stellar energy intercepted by the 


solids is Ld = Ad/^nal. If m m , n , m max , and q never change, the cross-sectional area Ad ,o and 


the initial dust luminosity Ld ,o are simple functions of N„ 
distribution. Thus, 

Ad, o 


and 


A d {t) = 


L d (t ) = 


1 + t/t c 
Ld, o 


and the parameters of the size 


( 2 ) 


. . (3) 

1 + t/t c 

At early times (t t c ), the cross-sectional area, dust luminosity and total mass in the disk 
are roughly constant. At late times (t '^> t r ), the area, luminosity, and mass decline as tr 1 
f Wvatt fc Dent 2002; Dominik fc Decin 2003). 


2.2. Destructive Collisions 

The simple relations in eqs. (l)-(3) hinge on maintaining a power-law size distri¬ 
bution with an invariant slope for particles with m m ,; n < m < Tn max ■ This outcome 
requires destructive collisions among equal mass objects. Collision outcomes depend on 
the ratio Q c /Qd, where Q* D is the collision energy per unit mass needed to eject half 


lision energy 

per unit mass (see also 

Wetherill & Stewart 

1993; 

Williams & Wetherill 

1994; 

Tanaka et al. 

1996; S 

tern & Colwell 

1997; Kenvon & Luu 

1999; 

O’Brien & Greenberg 

2003; 

Kobavashi & Tanaka 

2010). For impact velocity v. Q c = /iv 2 /2 (mi + m 2 ), where /j = 


1 Formally, collisions destroy two identical particles on the time scale 2 t c ; t c is then the time scale to 
destroy a single particle. 
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?77.i77i2/ {jri\ + m 2 ) is the reduced mass for a pair of colliding planetesimals with masses rri\ 
and m 2 . For equal mass objects, Q c = v 2 /8. 

Following standard practice, 


Qd = Qb r c b + Q a Pp r c 9 


( 4 ) 


where is the bulk component of the binding energy, Q g p g rc 3 is the gravity compo- 

nent of the bind i ng energy, and ry is the radius of a merged pair of p lanetesimals (e.g., 
Benz fe Asphaudll999t iLeinhardt et al.l 120081: iLeinhardt fc Stewartl 120091 1 . For icy objects, 


we adopt parameters - Qb 


10 5 erg g 1 


cm 


/3t , fj h as —0.40, Q g ~ 0.11 erg g 2 cm 3 , 


and j3 g 


1.3 


Burchell et al 


2005). 


(see also 

Davis et al. 

1985: 

Holsai 

jph 

1 7 | 

1994: 

Love & Ahrens 

Rvan et al. 

1999): 

Arakawa et al. 

2002: 

Giblin et al. 

2004: 


Setting Q c ps Q* d establishes constraints on the particles destroyed by an adopted colli¬ 
sion velocity. Among large particles with r > 1 km, the collision energy must overcome the 
gravitational component of the binding energy: r max = (v 2 /8Q g p)° “. When two particles 
with r < 1 cm collide, the impact kinetic energy must exceed the strength component of 
the binding energy: r min = (8 Qb/v 2 ) 0A . In between these two limits, Q* D is always smaller 
than Q c . Thus, collisions with velocity v destroy particles with r min < r < r max . Because 
the debris from these collisions produces an equilibrium size distribution between r mm and 
r max , continued destructive collisions maintain this size distribution. 


To relate these constraints to the properties of the planet and the central star, it is 
convenient to use the Hill radius, 


R h — 


Mp 

3M* 


1/3 




( 5 ) 


which establishes a volume where the gravity from the planet overcomes the gravity from the 
central star. For particles orbiting with random inclination in a spherical shell surrounding 
the pla net, the collision velocity is a simple function of the orbital velocity i>k, v = f v Vj< 
( Kennedy fe Wvatt 2011 b For a shell with a = rjx Rh, v = 3 1 / 6 (G /pi) 1 / 2 f v Mp^ 3 M^ 6 a p 1 ^ 2 . 


For Fomalhaut b, we adopt a set of fiducial param eters to e val uate t he collision velocity 


and other aspects of the collisional cascade (see also [Kenned y_& Wyatt 201ll). A 10 M ffi 
planet orbits a 1.9 M 0 star at a p ~ 120 AU. Satellites with p — 1 g cm -3 and r « 100 km lie in 


(e.g., 

Hamilton & Burns 

1992; 

Hamilton & Krivov 

1997; 

Toth 

1999; 

Slicn & Tremaine 

2008; 

Martin & Lubow 

2011 

): more compact configurations evolve too quickly. Clouds with (i) 
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typical ma ss x r i — MhJM p = 0.01 and (ii) collision velocities slightly larger than Keplerian, 
f v = 1.25 (Kennedy fe Wyatt 2011), then have orbital velocity 


v ~ 0.32 



/^x-1/2 / M p y /3 / Mi, \ 1/6 / a p \-i/2 
V0.27 V 10 M® / Vl.9M 0 ) V120 AU/ 


( 6 ) 


Substituting this velocity into our expressions for r min and r max '- 

1.54 „ -- / \ - 0.77 


100 


fv 


1.25 


P 


M n 


10 M® 


m_y 0 - 77 _ 

0.2/ ^1 g cm -3 

0 . 5 ! / \ 0.255 




1.9 M & 


120 AU 


0.11 erg g -2 cm 1 - 7 

- 0.77 


- 0.77 


km . 


(7) 


and 


1.8 x 10 


-4 


fv 


1.25 


-5 


v^y / 2 ( Q b 

0.2/ \10 5 ergg _1 


Mr, 


10 Mn 


- 5/3 


Mr 


1.9 Men 


- 5/6 


cm 

a 


, 0.4 


5/2 


120 AU’) 


5/2 


fim. 


( 8 ) 


Objects with intermediate sizes - r min < r < r max - have smaller Q* D than particles with 
r = rmin or r = r max . Radiation pressure typically ejects particles with r < r bl where 
r b ^ r min■ Thus, all particles with r < r max are either destroyed or ejected. 


For simplicity, many analytic models adopt a Q* D which is independent of radius. Eqs. 
[ZHH] justify this assumption: any Q* D which ensures the destruction of objects with some 
maximum radius guarantees that collisions will also destroy all smaller objects. 


2.3. Size Distribution 


With r max known as a function of Q g , deriving A^o requires values for r min and q. In nu¬ 
merical simulations of coll isional cascades , the slope of the equilibrium power-law size distri- 


bution is o « 3.5-3.7 (e.g., 

Dohnanvi 

1969: 

Williams & Wetherill 

1994: 

O’Brien & Greenberg 

2003: 

Kobavashi & Tanaka 

2010 

). In models where Q* D is constant with particle radius, q ~ 


3.5. When the bulk strength component of Q* D declines with radius, q ~ 3.7. We adopt q = 
3.5. 


In an optically thin swarm of satellites, radiation pressure sets r min . For dust grains 
orbiting the central star, radiation pressure removes particles smaller than the ‘blowout’ ra¬ 
dius, r b ~ (3L*Q pT ./87rcGpM*), where Q pr is the radiation pressure coefficient which accounts 
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for absorption and scattering (e.g., Burns et ah 1979). Fomalhaut has M* = 1.9 M® and 

7 /irn for icy grains with Q pr = 1 and p — 1 g cm -3 . 


20 L ®; thus, r b 


When particles orbit a planet, ejection depends on the orbital velocity of a parti¬ 
cle _around t he p lanet relative to the orbital velocity of the planet around the star (e.g., 
Burns et al 19791) . Defining f3 = F r /F g as the ratio of the radiative force to the gravita¬ 
tional force, radiation ejects particles orbiting a planet when (3 > /3q(v/v p ), where v is the 
orbital velocity of a particle around the planet, v p is the orbital velocity of the planet around 
the star, and /3 0 ~ 1/3 to 1. Thus, r b < {?>L+Q pr /%'KcGpM i <)(v p /f3 0 v). In physical units, 


n 


-(f) 



(9) 


More massive planets hold onto smaller particles. For particles orbiting at a fixed fraction 
of the Hill radius, r b is independent of a p . 

Setting r min = r b and integrating over a power-law size distribution with q = 3.5 yields 
the initial surface area. For convenience, we separate the linear dependence of Ad$ on the 
cloud mass Md into a linear dependence on XdM p \ 


Ad ,o — 1-4 x 10 


24 


Xd 

0.01 




r \ —V 2 

• m.n.cr. \ ' 


10 M® ) V100 knJ V 100 Fm 


- 1/2 


cm 


( 10 ) 


With these parameters, the initial surface area is roughly 10 times the observed surface area 
of a dust cloud in Fomalhaut b. 


2.4. Collision Time 


Deriving the long-term evolution of Ad requires a numerical estimate for the collision 
time. The simplest approaches adopt the lifetime of the largest particle against collisions 
with identical particles . More elaborate tre at ments include the impact of c ollisions with much 


smaller particles (e.g. 

, Wyatt et al. 

2007a.t 

>; Kobavashi & Tanaka 

2010; 

Kennedy & Wyatt 

2011: Kennedv et al. 

2011 

). Because the lifetime depends on a variety of relatively unknown 


parameters, we consider the time scale for destructive collisions among identical particles. 


To estimate the collision time, we again consider a spherical shell with semimajor 
axis a = ip Rr and thickness da = 772 a. The volume of this shell is V = 47 ra 2 < 5 a = 
47 r// 3 ? 72 a/T/p/( 3 M*). If the cloud consists of a monodisperse set of particles with mass m, 












































N = M ( i/rn. To express the mass of the cloud relative to the mass of the planet, we define 
Xd = M d /M p and use m = Xk pr^^/3. Then, N = 3x d M p jMpr'l nax . For mono-disperse par¬ 
ticles, the collisional cross-section is o = 47Tr^ aa ,. Thus, No = 3x d M p /2pr max . Combining 
all of these relations and defining t c = V/(2N 0 ov), 


t r = 


o 7/2 7/2 

27r r){ 7/2 p r max a p ' 
313/6 G i/2 fv Xd M l le M p /3 


( 11 ) 


When the cascade contains substantial mass i n particles with size s smaller than r mnx , the 


collision time is different from the t r in eq.fTThe.g., 1 

Wvatt et al. 

2007a 

b; 

Kobavashi & Tanaka 

2010; 

Kennedv & Wvatt 

2011; 

Kennedv et al. 

201 

!). For the fragmentation parameters used 


in eq. [4] and the collision velocity from eq. [6j collisions between one object with r & r max 
and another object with r > 0.1 r max destroy both objects. Cratering collisions with a much 
smaller object, r <C r max , eject roughly 2.7 times the mass of the smaller object. Thus, ev¬ 
ery collision removes mass from the largest object. Accounting for cratering collisions and a 
broader range of catastrophic collisions (i) allows larger objects to participate in the cascade, 
increasing r maX) and (ii) increases the rate large objects lose mass, shortening t c . 


Deriving the impact of these additional collisions requires integrating the collision rate 
over the size distribution. However, the range of sizes included in the integration depends 
on the collision rate. Although it is possible to construct an iterative solution, most inves¬ 
tigators simply set r max as a free parameter and derive the collision time for an r min set by 
the blowout radius ry Th e revised collision time is then ~ 0.01-4 t c flWvatt et al.ll2007al: 


Kobavashi fe Tanaka!12010c Kennedy &; Wvattll201l[) . However, this factor depends on v, 
Q* D , r mini and the details of the size distribution. For simplicity, we add a multiplicative 
parameter a (< 4) to our expression for the collision time. In §4.51 comparisons between this 


analytic model and our numerical results allow us to infer a for irregular satellite systems 
in Fomalhaut b. 


Converting the parameters in eq. [TT] to physical units, the collision time is: 


t c = 4.3 x 10 6 


(Vi\ 7/2 (V2 \ (Ji 


0.2 J 


0.5 

P 


1.25 


-1 


x d 

0.01 


-1 




1 g cm 


-3 


100 km 7 V120 AU7 


10 M® 

\ 7/2 


-1/3 


ilA 


yr 


1.9 M q 

( 12 ) 


More massive clouds, planets, and central stars shorten the collision time. Clouds consisting 
of larger, denser satellites orbiting planets with larger semimajor axes lengthen the collision 
time. Collisional cascades remove ~ 90% of the initial mass in ~ 10 collision times. Thus, 
the lifetime of the cascade is a significant fraction of the lifetime of an A-type star like 
Fomalhaut. 
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2.5. Example 


To illustrate the analytic model, we consider a simple example. With standard param¬ 
eters, M* = 1.9 M 0 , L* = 20 Lq, a p = 120 AU, rj i = 0.2, r/ 2 = 0.5, p = 1 g cm -3 , Q?, = 


IQ- 5 er g g- 1 cm 0,4 , /,, = 1.25, and r ma3 , = 100 km, eqs. ((21) and (flUl) yield the time evolution 
of the surface area. To compare with observations of Fomalhaut b, we adopt nominal values 
and uncertainties for the surface are a, A& ~ ljl a5 x 10 23 cm 2 and age, tb ~ 2001 2 qq Myr ( e -g-> 


MamaiekH2012t iKenvon et ah][2014] , and references therein). The relative surface area of a 


model satellite swarm is then Ad(t)/Ab. For satellite swarms with Xd = 0.01 and a range of 
masses for the central planet, Figure [1] compares the time evolution of the relative surface 
area evolves with observations of Fomalhaut b. 


For the adopted parameters, satellite swarms around super-Earths with M p = 30- 
100 M 0 match the data. Models with M p = 10 M 0 and M p = 300 M 0 almost match 
the data. Satellite evolution around lower mass or more massive planets do not match the 
data. 

Adopting other parameters leads to similar conclusions. Factor of two (ten) changes in 
v m ax (x<i) modify the relative surface area by ~ 25% at ages of 100-400 Myr. A modest, 
20% increase in rp increases t c by a factor of two at the expense of a 5% reduction in A di0 . 
This change yields a better match to the data for models with M p = 10 M 0 , at the cost of 
worse matches for models with M p = 100 M 0 . 


2.6. Advantages and Limitations 


The analytic model has several clear advantages. It is conceptually simple, easy to 
modify, and straightforward to calculate. The observable quantities A d and L d have obvious 
relationships to the physical parameters. Generating ensembles of debris clouds for plausi¬ 
ble variations in the physical parameters allows robust comparisons with large sets of data. 
Comparisons between data and models yield important insights into the evolution of cir- 


cumstellar debris disks and swarms of circumplanetarv satellite svstems (e.g.. 

Wvatt 

2008; 

Kennedv & Wvatt 

2010, 

2011; 

Kennedv et al. 

2011; 

Matthews et al. 

2014 

, and references 


therein). 


Despite its broad success, the model does not address several interesting issues in the 
time evolution of irregular satellites. In current theory, p lanets capture material from a cir- 
cums t ellar disk to supply t he circumplanetary sw arm (e.g., Koch fe Hansenlhoil ; Quillen et al. 


20121; iGaspar et al.l 120131: iNesvornv et al.l 120141) . If some captured objects have r > r r , 


these objects may accrete material from the swarm and reduce the lifetime of smaller particles 
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considerably. Dynamical interactions among growing large objects conld lead to significant 
scattering of particles within the swarm and interactions with larger satellites closer to the 
planet. 


The analytic model also does not allow time variations in r min , r max , and the slope 
of the power law size distribution. In a real collisional cascade, small particles gradually 
chip away at the larger objects. Significant reductions in r max shorten the lifetime of the 
cascade. Among smaller particles with r « r min , the typical collision time of 10 2 — 10 4 
yr (for Ad « 10 25 — 10 23 cm 2 ) is not much longer than the typical orbital period of ~ 
100 yr. Radiation press ure typically takes many orbital periods to eject small particles (e.g., 
Poppe fe Horanvi] 20 111 and references therein). Thus, r m j n might be significantly smaller 


than r;, at early times, enabling a larger initial surface area which declines more rapidly with 
time. In between r mm and r max , it takes many collision times to establish a power law size 
distribution with q = 3.5. If the initial size distribution is far from equilibrium, the early 
evolution of Ad might differ significantly from predictions of the analytic model. 


Addressing these issues requires numerical simulations. For a satellite swarm where the 
orbital elements are fixed in time, it is straightforward to conduct a suite of coagulation 
calculations to learn how time variations in r min , r max , q, and other physical parameters 
impact the long term evolution of the satellite swarm. We describe our numerical approach 
in §|3] and then discuss the results of the simulations in §|4j 


NUMERICAL MODEL 


To perform numerical calculations of the collisional evolution of an irregular satellite 
system, we use Orchestra , an ensemble of computer codes for the formation and evolution 
of planetary systems. Orchestra includes a multiannulus coagulation code which derives the 
time evolut ion of a swarm of solid objects orbiting a central mass flKenvon fc Brom levl 12004 
2008 . 2012). Although this code was originally designed to follow solids within a circumstellar 
disk, it is straightforward to modify the algorithms to track solids orbiting within a spherical 
shell. Orchestra also includes an n-body code which fo l lows the tra jectories and dynamical 


interactions of large objects (IBromlev fe Kenvonl 120061.1201 ll. 120131 1 . In these calculations, 


we disable dynamical interactions between coagulation particles and the n-bodies mediated 
by tracer particles. 
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3.1. Numerical Grid 


We conduct coagulation calculations of particles orbiting with semimajor axis a inside a 
spherical shell of width 5a around a planet with ma ss M p . Within this shell, there are M mass 
batch es with characteristic mass m*. and radius r*. flWetherill &; Stewartlll993l: iKenvon fc Luu 
19981 ). Batches are logarithmically spaced in mass, with mass ratio 5 = rrik+i/mk- Each mass 
batch contains W particles with total mass Mk and average mass rhk = M^/Nk . Particle 
numbers W < 10 15 are always integers. Throughout the calculation, the average mass is 
used to calculate the average physical radius fk, collision cross-section, collision energy, and 
other necessary physical variable s. As mass is added and removed from each batch, the 
average mass changes (iWethe r ill fe Stewart 1993). 


Numerical calculations with <5 > 1 lag the result of an ideal calculation with infinite 
mass resolution (see the Appendix). Simulations with 5 = 1.05-1.19 yield somewhat bet¬ 
ter solutions to the evolution of 10-100 km objects than calculations with 5 = 1.41-2.00. 
However, the evolution of the cross-sectional area of a swarm of solids is fairly independent 
of 5. To track the evolution of the size distribution reasonably well, we consider a suite of 
calculations with 5 — 1.19 (= 2 1 / 4 ). 


In these calculations, we follow particles with sizes ranging from a minimum size r min 
to the maximum size r max . The algorithm for assigning material to the mass bins extends 
the maximum size as needed to accommodate the largest particles. When collisions produce 
objects with radii r < r m j„, this material is lost to the grid. 


When the average mass in a bin exceeds a pre-set promotion mass m pro , the code creates 
a set of n-bodies with masses equal to the average mass in the bin. In these calculations, 
m pro = 10 24 g ( r max fa 600 km). Promoted objects are assigned a random semi-major axis, 
a pro , in the range (a — 5a, a + 5a), a random inclination sin i, a random orbital phase, and 
orbital eccentricity e = 0. For this first exploration of the evolution, we include the gravity 
of the central planet but ignore gravitational forces of nearby stars. 


3.2. Initial Conditions 

All calculations begin with a swarm of planetesimals with initial maximum size 0 and 
mass density p = 1 g cm -3 . These particles have initial number density n$ and total mass 
Mq. For the simulations in this paper, we consider two different initial size distributions for 
the planetesimals. To follow the analytic model as closely as possible, one set of calculations 
begins with a power law size distribution, n{r) oc r~ q and q = 3.5. To study whether our 
calculations produce this equilibrium size distribution, we begin a second set of calculations 
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with a mono-disperse set of planetesimals. 


3.3. Evolution 


The mass distribution of the planetesimals evolves in time due to inelastic collisions. 
All planetesimals have the same collision velocity, which is fixed at the start of each cal¬ 
culation. As summarized in Kenyon & Bromley (12004 . 2008). we solve a coupled set of 
coagulation equations which treats the outcomes of mutual collisions between particles in 
every mass bin. We adopt the particle-in-a-box algorithm, where the physical collision rate 
is navfg, n is the number density of objects, a is the geometric cr oss-sect ion, v is the relative 
velocity from eq. El and f g is the gravitational focusing factor (IWetherill fe Stewartl Il993l: 
Kenvon fe Luulll998l) . The collision algorithm treats collisions in the dispersion regime - 
where relative velocities are large - and in the shear regime - where relative velocities are 
small (IKenvon fe LuulIl998l: iKenvon fe Bromlevl 120141 1 . 


For swarms of satellites, all collisions are in the di spers ion regime, where we adopt a 
variant of the piecewise analytic approximation of Spaute et ah (119911 . see also Kenyon & 
Luu 1998; Kenyon & Bromley 2012). When collisions involve particles with r < 300 km, 
f g < 3-4. As satellites reach sizes of 1000_-2000 km, f g < 50. Compared to simulations 
where f g > 10 3 — 10 4 (e.g., Kenvon & Bromley 2008), gravitational focusing has a modest 
impact on the evolution. 

For each pair of colliding planetesimals and the collision energies Q c and Q* D defined in 
the mass of the merged planetesimal is 


m = mi + m 2 — m esc , 


(13) 


where the mass of debris ejected in a collision is 


Wiese 


0.5 (mi + m 2 ) 



(14) 


The exponent bd is a constant of order unitv (e.g.. 

Davis et ah 

1985; 

Wetherill & Stewart 

1993; 

Kenvon & Luu 

1999; 

Benz & Asphaug 1999; 

O’Brien & Greenbert 

j 200.4 

Kobavashi & Tanaka 

2010; 

Leinhardt & Stewart 

2012). Here, we consider bd = 1. 


To place the debris in the grid of mass bins, we set the mass of the largest collision 
fragment as 

Qc 


Q 


-&L 


m P 


Wlmax,d ^TlL, 0 


(15) 
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To explore the sensitivity of the evolution to this algorithm, we set mr ,,n = 0-2 and bi = 0 or 
1 (IKenvon fc Bromlevl 120081: iKobavashi &: Tanakall2010l; IWeidenschillind l2010l) . Lower mass 
objects have a differential size distribution n{r ) oc r~ q . After placing a single object with 
mass m maXt d in an appropriate bin, we place material in successively smaller mass bins until 
(i) the mass is exhausted or (ii) mass is pl aced in the smallest mass bin. Any material left 
over is removed from the grid (see also Kenyon & Bromley 2015 ). 


4. CALCULATIONS 

To examine the long-term evolution of satellite swarms, we consider a baseline model 
where the central star has mass 1.90 M® and luminosity 20 L®. Planets with M p = 1, 3, 
10, 30, 100, or 300 M® orbit with a semimajor axis a p = 120 AU. The spherical cloud has 
Xd = 0.01, 7]i = 0.2, and ip = 0.5. Particles within the shell have p — 1 g enr 3 , r min = 
100 pm and r max = 50, 100, 200, or 400 km. The initial size distribution is a power law 
with n{r ) oc r~ q and q = 3.5. To establish collision outcomes, we adopt the fragmentation 
parameters - Qb, Q g , f3b , and j3 g - summarized in 1 12.21 For the largest object in the debris 
we set = 0. 

To check the sensitivity of the calculations to these parameters, we vary one parameter 
and hold others fixed. In turn, we derive results for Xd = 0.001 or 0.0001; r min = 10 pm or 
1 mm; ( Qb,Q g ) = (5 x 10 4 , 0.055) or (2.5 x 10 4 , 0.0275); and b L = 1. Although we perform 
these calculations for all combinations of M p and r max in the baseline model, we focus our 
discussion primarily on calculations with M p = 10 M®. 

In describing the results of the simulations, we consider the long-term evolution of the 
size of the largest object ( 1 14.II) . the size distribution 11 14.21) . and the cross-sectional area 
(IJOJ). In many of the simulations, the masses of 2-6 objects reach the promotion mass. 
After promotion into the n-body portion of Orchestra, continued growth leads to strong 
dynamical interactions among a few n-bodies. Several examples of the long-term actions of 
the n-bodies illustrate their likely impact on the rest of the swarm and satellites orbiting 
closer to the planet (1 14.41) . In 1 14.51 comparisons with predictions of the analytic model allow 
us to clarify how the behavior of the swarm differs in the two approaches. 
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4.1. Evolution of the Size of the Largest Object 

4-1.1. Baseline model 

In the baseline model, the radius of the largest object r max depends on the mass of the 
central planet and the initial r max (Figs. [2H3]). When M p = 30-300 M®, the orbital velocity 
at 0.2 Rh is sufficient to destroy 100 km objects. Occasional catastrophic collisions among 
pairs of 100 km objects destroy them completely; the population of these objects gradually 
diminishes with time. More frequent cratering collisions with much smaller objects chip away 
at the mass in every large object. The average radius of these objects gradually declines 
from 100 km at t — 0 to 60-70 km at t — 1 Gyr. 

When M p = 10 M®, collisions almost shatter pairs of 100 km objects. Barring collisions 
with smaller objects, these objects grow slowly. However, cratering collisions with other 
satellites reduce their mass faster than collisions with the largest objects increase their mass. 
Thus, r max gradually declines with time. 

When M p = 1-3 M®, the orbital velocity is not sufficient to destroy 100 km objects 
(Fig. [21 lower orange and magenta curves). Cratering collisions are also insufficient to reduce 
their mass. Although collisions destroy all smaller objects, these objects grow slowly with 
time. After 1 Gyr, satellites orbiting a 3 M® planet reach radii of 400 km; satellites orbiting 
a 1M® planet have r max > 1000 km. In some cases, these calculations yield several n-bodies 
which interact dynamically. We discuss several of these outcomes in 1 )4.41 

In calculations with smaller r max , it is easier to destroy the largest planetesimals. Cra¬ 
tering collisions are also more efficient. For M p = 3-300 M®, r max declines more rapidly 
with time (Fig. [3]). By the end of the calculation at 1 Gyr, the largest satellites have r max 
= 2.5 km for M p = 300 M® up to r max = 35 km for M p = 3 M®. When M p = 1 M®, the 
collision energy is still insufficient to destroy 50 km objects. Thus, these objects grow slowly 
to ~ 300 km over 1 Gyr. 

When r max is initially larger than 100 km, it is easier for large objects to grow over time 
(Figs. [2H2])- Satellites with initial radii of 200 km (400 km) reach radii of 1000-2000 km 
in 0.3-1 Gyr around 1-10 M® (1-30 M®) planets. In several calculations, pairs of large 
satellites promoted into the n-body code scatter one another out of the grid. Scattering 
leaves behind a few much smaller objects, reducing r max within the grid. 

Collisional evolution with larger satellites orbiting more massive planets always reduces 
the size of the largest objects. For 300 M® planets, catastrophic and cratering collisions 
diminish the sizes of the largest satellites by 25% to 40%. The reduction in size is smaller, 
10% to 25% for 200-400 km satellites orbiting 100 M® planets. 


r\j 
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4-1-2. Initial cloud mass 

Changing the initial mass of the cloud has an obvious impact on the long-term evolution 
of the largest objects. When the mass of the cloud is smaller, the collision time is longer. 
Evolution is correspondingly slower. Thus, the sizes of the largest objects remain closer to 
their initial values. 

Fig. [4] illustrates the impact of initial cloud mass on particle size for satellites orbiting 
300 M® planets. When the initial r max is small, differences in the evolution are obvious. 
When x ( i = 0.01, it takes only 30 Myr for r max to decline to 20 km. A factor of 10 reduction 
in the initial cloud mass increases this time scale to 300 Myr. Another factor of ten reduction 
increases the time sale to 3 Gyr. 

As we increase the initial r max , the initial cloud mass has a smaller and smaller impact 
on the overall evolution. Although reducing the cloud mass increases collision times, large 
objects already have long collision times. Slow evolution simply becomes slower. Larger 
objects are also more impervious to collisional destruction; reducing Xd makes them even 
more impervious. 


4-1.3. Size of the smallest particles 

Not surprisingly, modifying the initial size of the smallest objects in the grid has little 
impact on the evolution of the largest objects in the grid (Fig. El). When M p = 10 M® and 
the initial r max is 50-100 km, the collisional cascade effectively destroys the largest objects 
in the grid. Collisions with the smallest particles in the grid remove little mass from the 
largest objects. Changing r min by a factor of ten has little impact on r max . 

When the initial r max is larger, the largest objects grow slowly with time. Large objects 
accrete little mass from the ensemble of objects with radii close to r m ; n . Changing r m ; n 
barely modifies the accretion rate. 

At late stages in the growth of large satellites, stochastic variations in the collision rate 
among large satellites produces large changes in r max . As large objects grow, the rate of col¬ 
lisions with other large objects declines and becomes more random. These random collisions 
produce large fluctuations in the time scale for objects to reach sizes of 1000 km (Fig. El 
top curves). When these objects are promoted into the n-body code, random dynamical 
interactions then yield random drops (or spikes) in r max . 
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4-1.4- Mass of the largest object in the debris 

How we distribute the debris into smaller mass bins also has modest impact on the 
evolution of the largest objects (Fig. [6]). In our baseline model with uil ,o = 0.2 and = 0, 
debris tends to fill bins with larger masses than in calculations with o = 0.2 and b^ — 1. 
When the cascade destroys large objects, the mass loss rate from the grid depends on the 
rate collisions transport mass to smaller and smaller objects. Thus, we expect calculations 
with bL — 1 to lose mass somewhat more rapidly than those with bL = 0. When larger 
objects grow, cratering collisions are less important; the exponent &£, then has negligible 
importance. 

Our results confirm these expectations. In the top curves of Fig. [6] 200 km and 400 km 
objects grow with time. Superimposed on a gradual rise in r max with time, stochastic vari¬ 
ations produce a few random increases in r max . Later, random dynamical interactions eject 
objects promoted into the n-body grid. The final radii of large objects is fairly independent 
of b L . 

In the lower curves of the figure, 50 km and 100 km satellites get smaller and smaller 
with time. Despite the somewhat larger mass loss of models with Ijl — 1, the final radii are 
independent of b L . 


4-1.5. Binding energy 

As discussed in (J2J Q* D - the binding energy of satellites - sets the size of the largest 
object destroyed in collisions with fixed impact velocity. Smaller Q* D allows collisions to 
destroy larger objects. 

Fig. [7] illustrates the impact of smaller Q* D on the evolution of 100 km and 400 km 
satellites orbiting 10 M® planets. When r ma x = 100 km, the Q* D in our baseline calculations 
is small compared to the collision energy. Collisions completely shatter these satellites. 
Although reducing Q* D makes them easier to destroy, the amount of debris lost in a collision 
is fairly similar. Thus, the evolution of 100 km satellites is independent of Q* D . 

In our baseline models, 400 km objects grow throughout the calculation. At late times, 
these objects reach sizes approaching 2000 km. Reducing Q* D has a clear impact on the 
growth of these objects. A factor of two reduction in Q* D prevents large objects from growing 
past 1000 km. Collisions among smaller objects are more destructive, removing objects from 
the grid more rapidly. With fewer objects to accrete, the growth of the largest objects stalls. 

Another factor of two reduction in Q* D completely halts the growth of 400 km objects. 
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Catastrophic and cratering collisions reduce these sizes of these objects by 35% to 40% in 
IGyr. 


4-1.6. Initial size distribution 

As a final example in this sequence, we consider the impact of the initial size distri¬ 
bution. When calculations start with a mono-disperse set of satellites, cratering collisions 
do not occur. If collisions between the large objects are catastrophic, the debris populates 
smaller size bins. Cratering collisions begin. Continued cratering and catastrophic collisions 
populate smaller and smaller size bins. Once all of the mass bins have debris, catastrophic 
and cratering collisions remove mass from all bins in the grid. Compared to our baseline 
calculations, these calculations have somewhat more mass in larger size bins. For the largest 
objects, collisions with larger small particles have a larger collision energy than collisions 
with smaller small particles. Thus, calculations with a mono-disperse set of satellites evolve 
somewhat faster than our baseline calculations. 

If large object collisions promote growth, there is less debris in smaller mass bins. 
Cratering collisions are less effective in filling smaller mass bins. Collisions between objects 
in these smaller bins are less frequent; mass loss from the grid is smaller. As these calculations 
proceed, there is much more mass in large objects (which are stronger) than in small objects 
(which are weak). The largest objects then grow faster. 

Our simulations confirm these expectations (Fig. IHJ). When the initial r max is 50 km, 
collisions are destructive. The largest objects get smaller and smaller with time. At late 
times, the satellites in the mono-disperse calculations are somewhat smaller than those in 
calculations with an initial power-law size distribution. 

When the initial r max is 100 km, the evolution is very sensitive to the initial size dis¬ 
tribution. In §[21 the analytic model suggests r max = 100 km for satellites orbiting a 10 M e 
planet. For this initial size, cratering collisions are critical. When they are present, r max 
declines with time. When they are not present, r max grows with time. The figure shows this 
dichotomy clearly. With no initial size distribution, 100 km objects grow slowly with time. 
With the power law initial size distribution, r max gets smaller and smaller with time. 

Among larger satellites with r max = 200 km or 400 km, the evolution also depends on 
the initial size distribution. With no cratering collisions among smaller objects in the grid, 
the cascade in fairly inactive. The mass in the grid is nearly constant in time. With more 
mass in the grid, the largest objects grow more rapidly. Promotion into the n-body grid 
occurs earlier; dynamical interactions are more severe. In these calculations, scattering of 
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n-bodies produced from a mono-disperse set of satellites leaves behind a few small objects 
which have suffered few collisions and are close to their original sizes. 


4.2. Evolution of the Size Distribution 


stant mass flow from the largest objects to the smallest objects (e.g.. 

Dohnanvi 

1969; 

Williams & Wetherill 

1994; 

O’Brien & Greenberg 

2003; 

Kobavashi & Tanaka 

2010 

). When 


r min ~ 0, catastrophic and cratering collisions maintain a power law cumulative size dis¬ 
tribution n(> r ) oc r~ qc with q c & 2.5-2.7. When r min > 0, cratering collisions between 
objects with r < r min and those with r > r m i n do not occur. Fewer collisions reduces the 
mass flow rate for r < r min . Mass then builds up in bins with r > r min . This excess of 
particles increases the rate of cratering collisions among larger particles, creating a deficit 
among these particles. Together, the excess and th e deficit produce a ‘wave’ in t he power law 
size distribution (e.g.. ICampo Bagatin et al.lIl994t lO’Brien fe Greenberd 120031 1 . Over time, 
continued collisional evolution tends to produce other waves among particles with larger 
radii. 


Addressing wave production in a numerical simulation requires an artificial extension 
of the size distribution to much smaller sizes. For an adopted size distribution for r < r min , 
it is possible to calculate the collision rate a nalytically an d to correct the mass flow rate 
for particles with r « r m i n (e.g., lO’Brien & Greenbe rg 2003]). For Fomalhaut b, however, 
radiation pressure removes particles with r < r mm on time scales much smaller than the 
collisional time. While some of these particles might lie on orbits which occasionally bring 
them back into the satellite swarm, most never return. Thus, real satellite swarms likely 
have wavy size distributions. 


In the rest of this section, we examine how ‘equilibrium’ wavy size distributions depend 
on model parameters. To discuss the time evolution of these size distributions, we derive the 
relative cumulative size distribution. At each r*. in the grid, the cumulative size distribution 
n(> r) is the number of objects with radius larger than r. To isolate the waviness about a 
power law, we define the relative cumulative size distribution 

n c ,rei = n{> r)/n 0 r~ qn , (16) 

where n 0 is a normalization factor. For these calculations, we adopt q n = 2 and normalize 
the relative cumulative size distribution to 1 at 10 km or at 100 km. 


By normalizing every relative cumulative size distribution at 10 km or 100 km, we 
suppress the natural evolution of rt 0 with time. In all calculations, n 0 follows the standard 
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evolution of the total mass in satellites: roughly constant at early times and then declining 
linearly with time at later times. In this section, we focus on the evolution of the shape of 
the size distribution. We return to the long-term evolution of the surface area in 1 )4.31 


4-2.1. Baseline model 


In the baseline model, the initial size distribution is a power law from r m ; n = 100 /mi 
to Tmax = 50, 100, 200, or 400 km. Large satellites tend to grow. Catastrophic and crater¬ 
ing collisions slowly reduce the sizes of satellites. For the adopted starting conditions, the 
collision time is 10-100 Myr. We expect the debris from destructive collisions to establish 
an equilibrium size distribution on this time scale. 


Fig. [9] illustrates the evolution of n Ct7 . e i for calculations with M p = 10 M 0 and r max 
= 100 km. This evolution is very different fr om standard predictions for a collisional cas¬ 


cade, oc r qr and q r ~ 0.5-0.75q (e.g.. IPohnanvil Il969t Williams fc Wetherilll Il994l: 


O’Brien & Greenberg 2003 : Kobavashi fc Tanaka 2Qloh . After only 0.1 Myr, n c>re i develops 
a characteristic shape, consisting of (i) a steep rise from r max to r « 50 km, (ii) a gradual 
rise from 50 km to 0.1 km, (iii) a steep drop from 0.1 km to 30 cm, and (iv) a steep rise 
from 30 cm to 100 /im. The rise from 50 km to 0.1 km has several distinct large-scale os¬ 
cillations; other fluctuations are small relative to the overall trend. After 1-10 Myr, n c ^ re i is 
independent of time except at large sizes r > 10-30 km. 


In this calculation, the lack of particles with r < 100 /mi produces the pronounced wave 
in the size distribution from 100 /mi to 0.1 km. With no very small particles (r < r min ), 
collisions remove objects with r & r min at a lower rate, producing the excess of 100-1000 /mi 
satellites. In turn, these objects remove larger particles at a faster rate, creating the large 
deficit at 1 cm to 10 m. At larger sizes, the waves are due to (i) stochastic collisions among 
the largest particles, which leads to an intermittent supply of debris among smaller particles 
and (ii) less frequent destructive collisions among 0.1-10 km particles and 1-100 m particles. 


The time scale to set up the equilibrium size distribution is many collision times for the 
smallest particles. In this example, the collision time for a typical small particle is vA^jV ~ 
10 4 yr. The main features of the wavy size distribution develop in just a few collision times, 
~ 0.1 Myr. By roughly 1 Myr, n c , re i establishes a distinctive pattern from 100 /mi to roughly 
10 km, which remains fixed for 1 Gyr. At the largest sizes, fluctuations in the collision rate 


2 In our convention, we have three power law slopes, q (differential power law), q c (cumulative power law) 
and q r (relative cumulative power law). These have a simple relationship: q ss q c + 1 ss q r + 3. 
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produce a time varying feature which slowly grows with time. 

Fig. [TQ] shows how the depth of the wave depends on the mass of the central planet. 
In this example, the minimum in n c , re i at 30 //m grows shallower and shifts to smaller sizes 
with decreasing planet mass. These changes are solely a function of the collision time. With 
t c oc Mp 7 (eq. [TT]), a factor of 300 in M p corresponds to a factor of 6.7 in the collision time. 
Thus, satellites orbiting a 300 M® planet experience roughly 7 times as many collisions 
over a fixed time interval as those orbiting a 1 M® planet. With fewer collisions, the size 
distribution of satellites orbiting 1 M® planets follows the initial power law more closely. 
Despite the longer collision time, collisions still establish the steep size distribution among 
the smallest particles and produce clear waves throughout the size distribution. 

For any mass of the central planet, the waviness within n C}1 . e i is also sensitive to the 
cloud mass. Gradually reducing the initial Xd produces the same progression in the depth 
and position of the deep minimum as in Fig. [T0J With t c oc xf 1 (eq. fTTj) . factor of 7 reductions 
in Xd for M p = 300 M® yield n Ctre i similar to the magenta curve in Fig. [10J With larger 
reductions, n Cjre ; more closely resembles the smooth power law of the initial size distribution. 

Changing the initial r max has little impact on n Cjre i at 0.1-1 Gyr (Fig. |TT]) . For r max 
= 50-400 km, the normalized size distributions at 100 Myr are essentially identical from 
100 /irn to 10-20 km. At larger sizes, the deviations from a power law depend on r max . 


4-2.2. Size of the smallest particles 

For calculations with fixed M p , Xd, and r max , changing r min has the same impact on 
n Ctre i as changing M p or Xd iFig. fl2|) . For the set of parameters in our calculations, collisions 
between a large particle with radius r; and a small particle with radius r s <C ri remove 2-3 
times the mass of the small particle from the large particle. When r m i n increases, collisions 
remove less mass from all remaining particles. Among 0.1 km and larger particles, the smaller 
mass loss has a fairly small impact on n C}1 . e i. However, less mass loss produces a steeper size 
distribution for particles with r r m i n and a larger deficit in particles at 10-100 cm. Larger 
r m i n also shifts the minimum in n CjT . e / to larger radii. 

When r min decreases, collisions remove somewhat more mass from all particles. The 
size distribution for the smaller particles becomes less steep and the pronounced minimum 
at 10-100 cm grows smaller. In our example, calculations with r m j n 10 /irn nearly eliminate 
the deep minimum in n CjT . e ; at 10-100 cm. This example retains the waves in n c ^ e i for r > 
1 cm and the steep gradient for r < 1 cm. 
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4-2.3. Mass of the largest object in the debris 

Distributing the mass in the debris differently also has a clear impact on n c ^ re i at small 
sizes. In our calculations, all collisions among particles with r < 1-10 km have Q c /Q*d > 1- 
Both objects shatter. When hi ~ 0, most of the debris is placed in bins with masses close to 
the mass of the original particles. When b^ ~ 1, more debris is distributed among particles 
with much lower mass. Spreading debris around more mass bins fills in the minimum in 
'ft'c^rel • 

Fig. [T3l illustrates this point. In calculations with 6^ = 0, there is a clear minimum in 
n Ctre i at 10-100 cm. When b L = 1, the minimum is less deep; the slope of the size distribution 
to smaller radii is shallower. At large radii (r > 0.1 km), 6^ has little impact on the n c , re i: 
all relative size distributions are wavy and roughly constant from 0.1-100 km. 

In both examples in the figure, the shape of the size distribution is independent of r max . 
All of the size distributions have small-scale fluctuations about the general trend, but these 
are small compared to the overall trends as a function of radius. 


4-2-4■ Binding energy 

Reducing Q* D has a similar impact on n C)re i as changing bi (Fig.lHj). Smaller Q* D leads to 
more ejected mass per collision. More ejected mass enhances the mass excess at the smallest 
sizes, steepening the size distribution. A larger population of smaller particles enhances the 
deficit at somewhat larger sizes. As with r min , changing Q* D changes the depth and the 
location of the deficit. Larger Q* D reduces the deficit and shifts it to smaller sizes. Smaller 
Q* d adds to the deficit and shifts it to larger sizes. 

Reducing Q* D also adds to the waviness of n c ^ re i at larger sizes. With more mass loss 
in every collision with much smaller particles, large particles distribute more mass among 
smaller mass bins. Truncating the size distribution at non-zero r m ; n creates a waviness in 
this mass loss, which is enhanced as Q* D is reduced. 


4-2.5. Initial size distribution 

To judge how the size distribution evolves when the calculations start from a mono- 
disperse set of satellites, we consider the baseline model with r min = 10 //m instead of 100 //rn. 
The smaller r m ; n give us a better view of the evolution of the population of small particles 
and yields a good comparison for models with and without an initial size distribution of 
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particles with radii smaller than r max - 

Fig.dSlilhistrates the evolution of a baseline model with M p = 10 M® and r max = 100 km. 
At the start of this calculation, collisions between pairs of 100 km objects produce satellites 
somewhat larger than 100 km and substantial debris. After 0.1 Myr, n Ctre i has a clear excess 
of particles with radii of a few km and some smaller particles. By 1 Myr, the excess has 
grown considerably; there is a substantial debris tail down to 1 m. In another 9 Myr, the 
debris populates the full range in allowed particles sizes and establishes a characteristic size 
distribution which remains fixed for the rest of the calculation. 

The equilibrium size distribution in Fig. [15] has the same features as in the baseline 
model. At small sizes (10 /irn to 1 cm), there is a steep and slightly wavy power law with 
slope q ~ 4.5 ( q c ~ 3.5; q r ~ 1.5). At larger sizes, the power law slope is closer to q ~ 2 with 
significant waves. Close to r maX) the slope again steepens. 

Fig. [T6] compares snapshots of the size distribution for calculations with (‘sd’) and 
without (‘no sd’) an initial power law size distribution among particles with radii smaller 
than r max . At 10 Myr, calculations with a mono-disperse set of large particles have more 
mass. In these calculations, the largest particles lose less mass through collisions with much 
smaller objects. Thus, they diminish in size less rapidly with time. For particles which are 
large enough to escape destruction, the extra mass in the largest particles allows them to 
grow more rapidly. 

At 1 Gyr, the n C)Te i for 10 /im to 30-50 km particles is independent of the initial size 
distribution. In the ‘no sd’ models, the largest objects reach sizes of 200-300 km in 1 Gyr. 
Satellites in the ‘sd’ models gradually lose mass and reach sizes of 70 km after 1 Gyr. 
Despite this difference, smaller debris particles produced in the collisions of the largest 
objects have an identical size distribution. Thus, the long-term evolution of the smallest 
objects is independent of the starting point. 


4.3. Evolution of the Cross-sectional Area 


For Fomalhaut b, observations cannot measure the size of the largest object or discern 
the size distribution across any range of sizes. Aside from the size and the color of the cloud 
as a whole, the only observable is the total brightn ess. For grains with radii > 10 /irn, we 


relate the brightn e ss to the cross-sectio nal area (e.g.. ICurrie et al.ll2013l: iGalicher et al.ll2013 


Kalas et al.l 120131: iKenvon et al.112014T) . To compare model results with the observations, 
we rely on the time variation of the cross-sectional area in each calculation relative to the 
adopted area for Fomalhaut b, A d 


10 23 cm 2 . A successful calculation matches this area at 
















the adopted age of Fomalhaut, tp ~ 200 Myr. We assign a factor of two uncertainty to A d 
and tp. The model ‘target’ is then a rectangular box in ( t F , A d ) space (see also Fig. [[]). 


4.3.1. Baseline model 

Fig. [T7] illustrates the time variation of A d for baseline models with initial r max = 
400 km. Initially, the surface area - A d oc M d oc x d M p - depends only on the initial cloud 
mass. Collisions then redistribute mass through the grid and gradually change the relative 
surface area. From eqs. dHanddU the collision time is t c oc Mp^Mj 1 . For fixed x d , swarms 
around more massive planets have larger cloud masses and shorter t c . Thus, calculations 
with M p = 100-300 M® evolve more rapidly than those with M p = 1-3 M®. After a brief 
re-adjustment where A d grows substantially, the relative surface area declines from ~ 100 
(where the swarm may become optically thick, see the discussion in 1 15.1.21 below) at 1 Myr 
to ~ 1 at 1 Gyr. Although the 100 M® models graze the target box at 400 Myr, these models 
generally fail to match the observations. 

Calculations with M p = 1-3 M® also fail. These models begin with relative surface 
areas close to the target and fairly long collisions times of 20-30 Myr. However, the largest 
objects in these calculations grow with time, removing small particles from the grid. After 
100 Myr of evolution, swarms of satellites orbiting 1-3 M® planets have relative surface areas 
at least a factor of two below the observations. 

Swarms of satellites around 10-30 M® planets match the observations throughout the 
100-400 Myr target period. In these systems, the initial cross-sectional area is roughly 10- 
20 times the area of Fomalhaut b. Although the largest objects in these simulations also 
grow with time, debris from the collisions of smaller objects maintains a large surface area 
for over 100 Myr. After 400 Myr (1 Gyr), satellites around the 10 M® (30 M®) planet have 
a surface area smaller than Fomalhaut b. Thus, there is a substantial cushion between the 
predicted time evolution and the evolution required to match the observations. 

Calculations starting with smaller r max yield smaller cross-sectional area at late times 
(Fig. fl8l) . With fixed total mass, swarms with smaller initial r max have larger surface area 
{A d oc r^Jx)- However, the collision time scales linearly with r max (eq.[H]). Thus, destructive 
collisions remove material more rapidly from swarms with smaller r max . More rapid mass 
loss results in smaller cross-sectional area. For M p = 30 M®, simulations with initial r max 
= 100-200 km still match the observations at 200 Myr. Satellites with initial r max = 50 km 
fall below the target. 

When r max < 200 km, calculations with M p = 10 M® also fail to match the observations. 
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Collisional evolution is too rapid for these models to match the observed surface area at 100- 
400 Myr. In contrast, models with more massive planets, M p = 100-300 M®, and r max = 
100-200 km, pass through the target. However, swarms with r max = 50 km succeed only for 
M p = 100 M®. 


4-3.2. Initial cloud mass 

The initial cloud mass has a more dramatic impact on the evolution of the surface area 
than r max . For fixed M p , the collision time scales inversely with Xd and r max . Smaller Xd 
and smaller r max lengthen the collision time. However, the initial surface area scales linearly 
with Xd and as rmax ■ Models with smaller Xd therefore start out with much smaller area and 
have more trouble matching observations. 

Fig. [EH] illustrates these points for M p = 10 M® and r max = 400 km. The baseline model 
with Xd = 0.01 matches the data for ages of 100-400 Myr. Factor of ten smaller swarms have 
an initial surface area somewhat larger than Fomalhaut b, but collisions reduce the area by 
almost a factor of ten after 100 Myr. Another factor of ten reduction in Xd leaves the surface 
area below observations throughout the evolution of the swarm. 

Fig. [2m shows that low mass swarms around massive planets can also match observations. 
When M p = 100 M®, satellites with Xd = 0.01 have a surface area that grazes the upper 
edge of the target box at 300-400 Myr. Reducing the initial mass by a factor of 100 yields 
a surface area that grazes the lower edge of the target box at 100-200 Myr. Intermediate 
cloud masses match the data well; models with Xd = 0.001 pass through the upper middle 
of the target box. 


4-3.3. Size of the smallest particles 

In the analytic model, the surface area of a satellite swarm scales with r nnn . Changing 
r m i n by an order of magnitude thus modifies the initial surface area by a factor of ~ 3. 
Because r min has little impact on the cloud mass or the collision time, the early evolution of 
satellites with different r min is nearly identical. Over time, however, the smallest particles 
shape the size distribution at larger sizes (Fig. [T2lh When r m i n is smaller (larger), more 
(fewer) intermediate particles with r fa 10 cm to 10 m survive the cascade. If these particles 
contribute much to the total surface area of the swarm, then we expect the surface area at 
late times to scale more steeply with the minimum particle size. 

Fig. I2T1 demonstrates that transformations to the size distribution have little impact on 
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the evolution of the surface area. In the baseline model with M p = 10 M®, r max = 200 km, 
and r min = 100 /mi, the evolution of the surface area passes through the lower left corner of 
the target area. In calculations with r min = 1 mm, the surface area is \/l0 smaller and falls 
well below the target. Factor of ten smaller r m ; n yields a factor of \/T0 larger surface area 

_-j/2 

which lies in the upper half of the target. Thus, the surface area scales exactly as r mi ' n . 

Despite significant differences in the size distributions of calculations with different r m , n , 
these changes have little impact on the evolution of the surface area. Particles with radii of 
1 cm to 100 m have a limited fraction of the mass of the larger particles and a negligible 
surface area compared to much smaller particles. Augmenting (or reducing) the population 
of these particles by factors of 100-1000 has no observable impact on the total surface area 
of the cloud. 


4-3-4■ Mass of the largest object in the debris 

In our collision model, we use a simple algorithm to distribute debris into the mass 
bins. For the baseline model with bi = 0, the largest object in the debris has 20% of the 
total mass in the debris. In comparison models with b^ — 1, larger relative collision energies 
place more material in lower mass bins. Material leaves the grid more rapidly. Although the 
distribution of debris has little impact on the evolution of the largest objects (Fig. [6]), the 
relative number of the smallest particles depends on the mass of the largest object in the 
debris (Fig. [13]). In the complete ensemble of calculations, N(r min ) changes by a factor of 
3-5; we expect similar changes in the cross-sectional area. 

Fig. ED confirms this conjecture for calculations with M p = 10-30 M®, Xd = 0.01, r max 
= 400 km, and r min = 100 pm. When b L = 0, the predicted surface area passes through the 
target box. For b L = 1, calculations with M p = 30 M® pass through the lower edge of the 
target; models with M p = 10 M® completely miss the target. 

Larger bi allows satellite swarms around more massive planets to match observations 
at 100-400 Myr. In the baseline model, calculations with M p = 100-300 M® and Xd = 0.01 
have much larger surface area than Fomalhaut b. When 5® = 1, satellites orbiting planets 
with M p = 100 M® (300 M®) pass through the lower (upper) half of the target box. 


4-3.5. Binding energy 

In our suite of calculations with different Q* D , smaller Q* D prevents the largest objects 
from growing (Fig. [7]). The additional debris produced from collisions between the largest 
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objects decreases the relative numbers of satellites with radii larger than 1-10 cm (Fig. [T5]) . 
However, smaller Q* D has little impact on the relative population of the smallest objects 
which dominate the cross-sectional area. Thus, swarms of satellites with different binding 
energies have comparable surface areas. 

Fig. 1221 shows the modest variations in cross-sectional area for calculations with different 
Q* d . In the baseline model, 400 km satellites reach radii of 2000 km on time scales of 100- 
400 Myr. Large fluctuations in the surface area resulting from stochastic variations in debris 
production from the collisions of the largest objects begin at 100 Myr and continue beyond 
1 Gyr. Lowering Q* D by a factor of two limits the growth of these objects to 1000 km. 
Slower growth delays the large oscillations in the surface area and minimizes them at later 
times. Another factor of two reduction in Q* D completely eliminates the growth of the 
largest objects. Less growth enables more debris and larger surface area at early times. As 
the calculation proceeds, debris production is somewhat more modest than calculations with 
larger Q* D) resulting in a slightly lower surface area at later times. 

In calculations with b L = 1, smaller Q* D speeds up the time evolution of the surface area. 
When Q* d is smaller and = 1, collisions place more debris in smaller mass bins. Placing 
debris in smaller mass bins allows the cascade to eject mass from the grid more rapidly. 
Thus, the surface area declines more rapidly with time. For the smallest Q* D considered in 
our calculations, models with 5^ = 1 require very large and probably unrealistic initial cloud 
masses to match the target for Fomalhaut b. 


4-3.6. Initial Size Distribution 

When calculations begin with an initial power-law size distribution, the initial surface 
area is substantial (Fig. [TTlh With a mono-disperse set of large particles, the initial surface 
area is negligible. For r max = 100 km and r min = 100 /mi, it takes several collision times to 
populate the low mass end of the size distribution (Fig. (THD and increase the surface area. 

Fig. [23] compares the growth of the surface area for swarms of satellites with (‘sd’) and 
without (‘no sd’) an initial power-law distribution of small particles. When r max = 100 km, 
the collision time is rather short. It takes only a few x 10 6 yr for the surface area to reach the 
levels of the baseline model. At 10 Myr, the satellite swarm in the mono-disperse model has 
relatively more small particles than the swarm with the initial power-law size distribution 
(Fig. fl6l) . Thus, the mono-disperse model has a larger surface area. By ~ 1 Gyr, collisions 
have erased the starting conditions; the size distributions (and hence the surface areas) are 
indistinguishable. 
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When the initial r max is 400 km, the evolution of satellite swarms with and without an 
initial power law size distribution diverge. In these calculations, the collision time is much 
longer. Populating the small end of the size distribution takes 30-40 Myr instead of 3-6 Myr. 
With fewer small objects to chip away at the mass of the largest objects, the largest objects 
grow more rapidly and contain a larger fraction of the initial mass (Fig. [8j). At ~ 100 Myr, 
these calculations produce more debris than the baseline model. As the evolution proceeds, 
the larger mass tied up in the largest objects speeds up the decline of the population of 
small particles. After 400-500 Myr, the surface area in the small particles is a factor of > 2 
smaller than in the baseline model. 


4.4. Growth of Very Large Objects in the Swarm 


In some calculations, the largest objects in the swarm grow by factors of 10-1000 in 
mass over 0.1-1 Gyr. As these objects grow, their gravity may stir up smaller objects in the 
swarm. Small amounts of stirring increase typical collision velocities and enhance the impact 
of the collisional cascade. Larger amounts can eject material from the swarm. Aside from 
ejections, dynamical interactions between pairs of very large objects can place satellites on 
bound orbits closer to the planet. If the planet already has a set of closely bound satellites, 
a ‘new’ satellite might disrupt the indigenous satellite system. 

To explore the impact of these processes within our simulations, we begin with general 
principles. When large satellites grow, they try to stir up smaller satellites to their escape 
velocity. Satellites with radii r max have escape velocity 

Vesc = 0-9 ( ) kms -1 . (17) 

VlO 3 km/ V ; 

For a central planet with M p = 10 M®, satellites with r max & 350-400 km have an escape 

velocity, 0.3 km s -1 , comparable to the collision velocity of satellites in a spherical swarm. 

Thus, stirring is a crucial issue for large satellites with radii of 1000-2000 km. 


Satellites orbiting the planet have Hill spheres where the gravity of the satellite domi¬ 
nates the gravity of the planet. For satellites orbiting at semimajor axis a, the Hill radius is 
Rh,s ~ r .s a (e.g., eq. [5]) where 


0.04 


1000 km/ \10 M® 


M n 


- 1/3 


(18) 


Larger satellites around less massive planets have larger Hill radii. 


When the semimajor axes of orbiting satellites differ by 3-4 Hill radii or less, they 
interact dynamically. For simplicity, we define the minimum orbital separation for stability 
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as A a « 4 r s a. In a system of N massive satellites within a spherical shell of width 5a = 772 a, 
the system is dynamically stable when 4 Nr s a < 772 a. Solving for N, the maximum number 
of non-interacting satellites is 



(19) 


When satellites orbit massive planets with M p > 1-300 M 0 , large-scale dynamical interac¬ 
tions require a few satellites with r max > 1000 km. For fixed r max , dynamical interactions 
are more common around less massive planets. 

In our suite of calculations, stirring of smaller satellites by the largest satellites is a 
relatively minor issue. For objects within the collisional cascade, the initial collision velocities 
produce shattering. Larger relative velocities produce a little more debris and a somewhat 
faster decline in the relative surface area of small particles. Larger objects are already 
immune to the cascade; larger relative collision velocities have little impact on their evolution. 

Dynamical interactions among 77 -bodies are more important. As one example, Fig. 1251 
tracks the time evolution of the semimajor axes for a set of 'n-bodies orbiting a 1 M 0 planet. 
In this calculation, the initial set of 400 km objects grows throughout the evolution of the 
satellite swarm (Fig. |2J). At ~ 60 Myr, the coagulation code promotes the first satellite into 
the n-body code with a circular, but highly inclined, orbit at 0.15 AU. Roughly 10 Myr later, 
a second ?r-body appears with an orbit at 0.18 AU. Within another 5 Myr, a third 77 -body 
has an orbit at 0.21 AU. At 100 Myr, the separations of these satellites are roughly 3.5 
mutual Hill radii. Strong dynamical interactions are inevitable. A scattering event between 
the outer two satellites places one on a very close orbit with the inner satellite. All develop 
eccentric orbits. Eventually, the more massive inner satellite ejects the other two satellites 
and ends up on an eccentric orbit much closer to the planet. 

In this suite of calculations, the dynamical evolution of the 77-bodies has no impact 
on satellites remaining in the coagulation code. A few larger objects continue to grow. 
Promotion of two of these satellites into the 77-body code leads to another set of dynamical 
interactions at 200 Myr, where the two new (and lower mass) 77-bodies are ejected and the 
original massive 77-body moves a little closer to the planet. One last satellite promoted into 
the 77-body code at ~ 500 Myr orbits on a circular, highly inclined orbit well away from the 
inner massive satellite. This system remains stable for the rest of the calculation. 

At the end of this calculation, roughly 33% of the initial mass in solids remains in orbit 
around the planet. Nearly all of this material is in the two large satellites orbiting at 0.05 AU 
and 0.19 AU. Dynamical interactions (44% ± 7%) and radiation pressure (56% ± 6%) eject 
equal amounts of material. Despite this rough equality in mass, radiation pressure removes 
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100 (.i m particles from the grid. Dynamical interactions place four Pluto-mass planets into 
orbits around the central star. 

This evolution of large objects is fairly typical. Roughly half of the simulations with 
growing n-bodies leave massive satellites orbiting the planet. Nearly 15% of these have 
satellites orbiting at semimajor axes well inside the initial extent of the satellite swarm. Less 
than 5% have satellites outside the initial boundary of the swarm. In the rest, 1-2 satellites 
orbit stably within the swarm. 


4.5. Comparing the Analytic and Numerical Models 

Compared to the analytic model, the numerical simulations yield several clear differ¬ 
ences in the collisional evolution of a satellite swarm. The sizes of the largest objects change 
considerably in 0.1-1 Gyr. When catastrophic and cratering collisions dominate, r max de¬ 
clines by 30% or more. The largest objects then have roughly 30% of the mass of the largest 
objects in an analytic model where r ma x is constant in time. 

In some systems, the largest satellites grow substantially. Left unchecked, this growth 
yields massive objects capable of disrupting the satellite swarm (and perhaps satellite systems 
closer to the planet). Swarms orbiting lower mass planets are more prone to this evolution 
than swarms around more massive planets. 

The size distribution does not follow a simple power law. Although the numerical 
simulations roughly follow this power law for satellites with r > 0.1 km, all of these models 
produce large (factor of 3-5) waves about the power law. At small sizes, there is a large 
deficit in 0.1 cm to 10-30 m particles relative to the analytic power-law. In calculations with 
an initial power-law size distribution, smaller r m ; n and larger b l reduce the size of the deficit. 
Calculations starting with a mono-disperse size distribution also yield smaller deficits. 

The smallest particles with r < 0.1 cm follow a steep power law with q & 4-5. In 
our simulations, the lack of particles with r < r min limits mass loss among larger particles. 
When these particles lose mass less rapidly, they stay in the grid for longer periods of time, 
steepening the size distribution. 

Despite these differences, the time evolution of the surface area in the baseline model 
(Fig. 1171) is similar to the time evolution of the analytic model (Fig. [[]). In the analytic 
model, satellite swarms orbiting 10 M® planets evolve through the middle of the target. 
Swarms around 30 M® planets have a surface area in the upper half of the target box. 
Swarms orbiting smaller (3 M®) or larger (30 M®) planets have areas that graze or just miss 
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the target. In the baseline numerical model, satellites orbiting 10-30 M® planets match the 
observations. Swarms around 3 M® or 100 M® planets fall below or graze the upper edge 
of the target. Overall, the numerical simulations require slightly more massive planets to 
match the observations than the analytical models. 

Fig. [26] compares the evolution of the mass in a baseline model with the predictions of 
the analytic model using three different values for the correction factor a. For M p = 10 M®, 
Xd = 0.01, r max = 100 km, and r min = 100 /mi, the mass in the numerical model begins to 
decline at 0.1 Myr. It takes 40-50 Myr for the mass to reach 10% of the initial mass and 
700-800 Myr to reach 1% of the initial mass. At late times, the mass evolves with time as 
M(t ) oc t~ n and n fa 0.8-0.9. Thus, the mass declines somewhat less rapidly than in the 
analytic model. 

In the analytic model, the mass loss rate is initially smaller than in the numerical model. 
As time proceeds, analytic models with a = 0.333 (3.333) decline more rapidly (slowly) than 
the numerical model. When a = 1, the analytic and numerical models match at 20-30 Myr. 
After this time, the t _1 decline of the analytic model results in a faster rate of mass loss 
than the numerical model. 

For comparisons between the analytic model and the complete suite of numerical sim¬ 
ulations, the ‘best’ a depends on the input parameters. For all models with the baseline 
set of parameters, a fa 1. Changing Xd, r min , Q* D , and the initial size distribution has little 
impact on a. When &£ fa 1, the faster removal of material in the grid leads to more rapid 
evolution and smaller a fa 1/2 to 1. 

5. DISCUSSION 

Our suite of simulations paints an interesting picture for the evolution of swarms of 
satellites orbiting 1-300 M® planets at 120 AU from a central 1.9 M 0 star. Depending on 
M p , Xd, and r max , the masses of the largest satellites grow or shrink with time. Satellites 
with large Q* D grow; satellites with small Q* D shrink. Clouds with smaller Xd and larger r max 
evolve more slowly. Outcomes are insensitive to r m j n , the shape of the initial size distribution, 
or the algorithm for distributing debris among the mass bins. 

However the largest objects evolve, the cumulative size distribution transforms into a 
standard shape which is fairly independent of the input parameters. This standard shape 
has three distinct pieces: (i) a steep power law with n(> r) oc r~ qc and q c fa 4 ( q r fa 2) 

for small particles (r < 10-30 cm), (ii) a flat portion where q c fa 0-1 ( q r fa —2 to —1) 

for intermediate size particles (r fa 30 cm to 0.1 km), and (iii) a shallow, wavy power law 
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with q c fa 1-2 (q r fa —1 to 0) for large particles (r > 0.1 km). The development of this 
standard shape depends on the collision time: swarms with longer collision times take longer 
to establish this size distribution. 

Despite the diverse outcomes, a broad set of satellite swarms produce a cross-sectional 
area which matches the observed area in Fomalhaut b. For models with the nominal r m j n , 
Q* d , and bj^ Fig. [271 summarizes these outcomes a function of M p , r max , and Xd- Swarms 
orbiting low mass planets always fail. Satellites around more massive planets are often 
successful. Successful models have a factor of 2-3 range in the initial mass in satellites 
relative to a ‘best’ model with an initial mass of rsj 10 27 g. 

Aside from the initial mass, the size of the smallest stable particle orbiting the planet 
and the mass distribution of debris from a collision establish the ability of a model to match 
the observed cross-sectional area. If particles with r min = 10 /irn can stably orbit the planet, 
smaller initial cloud masses are possible. Larger r min and > 0 reduce the grid of successful 
models. 


5.1. Theoretical Issues 

In these coagulation calculations, the standard size distribution has several features 
in common with results for debris disks orbiting main sequence stars. There are also a 
few major differences. In addition, several uncertain parameters establish whether model 
satellite swarms have cross-sectional areas at 100-400 Myr comparable to the observed area 
in Fomalhaut b. Assigning different values for r min , Q* D , and 6^ allows swarms with different 
combinations of M p , Xd, and r max to match the observations. Here, we discuss features of 
the size distribution and consider the flexibility of the theory in setting the various input 
parameters and the likely consequences of our choices. 


5.1.1. Size Distribution 


In Figs. 171117)1 the size distributions derived in our calculations are radically different 
from the smooth p ower law, n(r) oc r~ q with q sa 3.5-3.7, expected from an equilibrium 


collisional cascade flDohnanvilll9fi9t lO’Brien fe Greenberdl2003l: iKobavashi fe Tanakall2010l) . 
The general trend of n(r ) is much shallower than this power law. Pronounced waves are 
superposed on this general trend. 

Wavy size distributi ons are a common feature in numerical calculations of debris disks or- 
biting 1-3 M e stars (e.g., Campo Bagatin et all 1994 : O’Brien fe GreenbereJl2003 : Thebault et al 
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20031; iKrivov et al.l 120061: Ihdhne et al.l l2008t iGaspar et al] 12012c Krai et al.l 12013 ). These 
waves have two typical sources: the low mass cutoff of the size distribution and the transi¬ 
tion between the bulk strength and gravity regimes in analytic expressions for Q* D (see eq. QJ. 
Although the amplitudes (in n{r )) of the waves depend on the radius of the low mass cutoff, 
the ratio Q c /Q*d, and the parameters in the relation for Q * D , typical values are a factor of 
10 or smaller. 

By typical debris disk standards, the waves in Figs, [HI 17711 are somewhat extreme. For 
our calculations, the positions and relati ve spaci ng of minima and maxima follow general 


predictions from analytic models (e.g.. lO’Brien fe Greenberg 2003). Factor of 3-10 ampli¬ 
tudes at 0.1-100 km are also normal. However, the amplitudes of waves at small sizes are 
10-1000 times larger than those reported from numerical simulations of debris disks. 

The long-term collisional evolution of satellite swarms is responsible for this difference. 
Compared to debris disk s ar oun d stars, sa tellite swarms orbiting planets are much more 
collisionally evolved (e.g.. IBottke et al. 2 0101 ). More collisional evolution enhances the excess 
of particles with r ~ 1-10 r min and the deficit of particles with r « 10-100 r min . As a result, 
the amplitude of the wave simply grows larger and larger with time. Fig. [TU] clearly shows 
the impact of longer collisional evolution: massive systems with shorter collision times have 
much stronger waves than low mass systems with longer collision times. 


5.1.2. The minimum particle size 


In the baseline model, r min is the blowout radius r b for a 10 M® planet. The slow 
variation of r b with M p justifies this assumption for all M p (eq. [HI). However, radiation 
pressure cannot remove small particles when (i) the collision time is comparable to or shorter 
than the orbital period around the planet and (ii) the optical depth of the cloud is one or 
larger. If satellite swarms meet either of these conditions, smaller particles stably orbit the 
planet. 


The optical depth r of the cloud is roughly the ratio of Ad to the cross-sectional area 
defined by the physical extent of the swarm A s = Setting r = Ad/A s and adopting 

the nominal parameters, 


9 x 10' 4 ( —A: 
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( 20 ) 


For planets with M p < 30 M®, satellite swarms have Ad < 10 25 cm 2 throughout their 
evolution. These swarms are never optically thick. Early in the evolution of swarms orbiting 
more massive planets, Ad > 10 26 cm 2 ; these swarms are optically thick. 
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For a single small particle, it is straightforward to derive the ratio £ of the collision 
time to the orbital period. The collision time is roughly t c fa V/vA the orbital period is 
T = 2na/v. The ratio is then: 


£ fa 1.85 x 10 ( 
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( 21 ) 


Around low mass planets with M p < 10 M®, the initial cross-sectional area of the swarm 
is < 10 24 cm 2 . Radiation pressure removes small particles faster than collisions. Among 
more massive planets, A^ > 10 26 cm 2 for t < 1 Myr. High speed collisions destroy small 
particles faster than radiation pressure removes them. The smallest particles in the swarm 
are then much smaller than 100 //m. Although the cross-sectional area of these swarms is 
then formally very large, the observed area is limited by the optical depth. With r ~ 1, the 
maximum area is roughly 10 26 cm 2 . 

This discussion implies that r min is rarely smaller than the nominal blowout size ry. 
Early in the evolution of satellite swarms around massive planets, particles with sizes smaller 
than 10-100 /jm remain bound. This phase is short-lived, ~ 1 Myr. As these systems evolve, 
Ad declines rapidly as r min returns to its nominal value. The evolution then continues as 
outlined in 1 14.31 


5.1.3. The particle strength 


The binding energy of solid particles establishes collision outcomes. When the colli¬ 
sion energy Q c exceeds the binding energy Q* Dl more than half of the mass of the colliding 
pair of particles ends up in debris. In circumstcllar disks, collisional damping and gravi¬ 
tational interactions between particles often limit the impact of Q* D on the evolution (e.g., 
Kenvon fe Bromley 2008. 2010). For satellite swarms within a spherical shell, damping and 
gravitational interactions are minimal. With Q c solely a function of r/i and M p , the growth 
of the largest particles is a strong function of Q* D . 

For icy objects with sizes much larger than 1 cm, astronomical observations, laboratory 
experiments, and numerical simulations paint a disparate picture for Q* n . Recent experi- 
ment s colliding cm-sized icy solids in the lab sugge st tensile strengths of roughly 10 6 erg g _1 


fe.g.. |Shimaki feArakawal 120121; lYasui et al .1120141) . which agrees with previous results (e.g., 
Rvan et al.l Il999h . Numerical simulati ons of high speed collis i ons between icy objects ar e 
the basis for the expression in eq. 0] fjBenz fe Asphaud Il999l: iLeinhardt fc Stewartl 120091 1 . 
Typically, Qb m 10 5 — 10 8 erg cm 0 4 g _1 and Q g = 0.1-2 erg cm 1,65 g _1 . For r fa 1-10 cm, the 

Experiments and simulations 


simulations suggest a binding energy of roughly 10 6 erg g 1 
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agree on the strength for small objects. 


Observations of comets yield much smaller binding energies. Models for comet D3/1993 
F2 ( Shoemaker-Levy) and o t her disrupted comets suggest binding energie s of 1-10 3 erg g _1 


(e.g. 


Asphaug fe Benzlll996t iRichardson et al.ll2007l; ISkorov fe Bluml 12012 ). Data from Deep 


Impa ct imply a strength in the middle of this range (e.g.. lHolsapple fe Housenll2007l: I A ’Hearn 


2011). With r fa 0.1-1 km for the nuclei of these comets, the maximum strength of 10 3 erg g 


-i 


is a factor of 10 or more smaller than expected from eq. Q] and the results from laboratory 
and numerical experiments. 


Our adopted values for the parameters in Q* D lie intermediate between observations and 
numerical simulations. Values for Q* D similar to results from the studies of comets preclude 
the growth of large satellites around planets with M, p = 1-30 M® (see Fig. [2]). More material 
then participates in the collisional cascade. Although evolution times are somewhat shorter, 
the evolution of the cross-sectional area is unchanged. Thus, significantly smaller Q* D does 
not change our conclusions. 


Larger values for Q* D enha nce t he grow th of large satellites around all planets. For 
sufficiently large Q* D as in iBenz fe Asphaug! (19991), the collisional cascade is limited. Few 
satellite swarms have sufficient surface area to match observations of Fomalhaut b. In these 
systems, the largest satellites probably grow large enough to disrupt the satellite swarm 
completely (see Tl.4p . Then, all models fail: none match observations of Fomalhaut b. 


5.1.4- The size of the largest particle in the debris 


Theory currently provides limited guidance on m m ax,d, the mass of the largest par¬ 
ticles in a cloud of debris ejected during a high velocity collision. For cratering colli- 
sions- m Psr is a simple fun c tion of the co llision energy and the gravity of the planet (e.g., 
Housen fe Holsapple boil : Svetsoy 2011 ). However, there are no direct calculations of 


Wcthcrill 


fe StewartT f 1993 ) examined laboratory data and adopted the simple rela- 


Wlmax,d- 

tion / //„,„ ,■,/ = 0.2m f 9r used in our baseline mo del. Recent experiments confirm this choice 


(e.g. 


Burchell et al.ll2005t iPoelchau et al.ll20l4 and references therein). 


For catastrophic im pacts, numerical simulations provide somewhat conflicting advice 
for m maX 'd- hi Lei nha r dt fe Ste w a rt (2012), collisions of 1-10 km icy objects yield a broad 
range, rn max ,d/m esc ~ 0.001-1. Power-law fits to the distribution of debris particles require 
assigning either a lower bound to the size of a debris particle or a slope to the size distribution. 
Choosing the slope leads to a fixed value rn maXt d/rn esc ~ 0.026 independent of Q* D . 
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Durda et ah (2004, 2007) describe numerical simulations of collisions for 10-100 km 


he debris is 1-2 


rocky objects. In these simulations, the mass of the largest object within 
orders of magnitude sma ller and very sensitive to the ratio of Q c to Q* D (IMorbidelli et ah 


20091: iBottke et al.ll2010[). These ca lculations predict much steeper size distributions than 


those of Lei nhard t fe Stewartl (2012). However, including these results in coagulation codes 
requires adopting a shallow slope for small sizes to conserve mass which introduces additional 
input parameters. 

For high velocity collisions of small objects, laboratory experiments suggest the mass 
in the largest debris particl e scales roughly inversely with the ratio of the collision energ y 
to the binding energy (e.g., lArakawalll999t lArakawa et al.ll2002l; IShimaki fe Arakawall2012l) . 
Available data imply larger particle sizes for collisions between more porous and stronger 
targets. 


Our approach to placing debris in mass bins roughly follows the spirit oflLeinh ardt & Stewar t 


(2012). In our baseline model, debris from cratering collisions agrees with experimental re¬ 


sults; the mass of the largest object in catastrophic collisions is roughly 10 (ILeinhardt fe Stewart 


2012 ) to 100 (IBottke et ah 20101) times larger than inferred from numerical simulations. 
Compared to the predictions of numerical collision codes, these calculations probably under¬ 
estimate the rate of decline for the cross-sectional area around massive pla nets. In models 


with = 1, catastrophic collisions yield results more similar to Leinhardt & Stewart (2012). 
Although this treatment of cratering collisions leaves too much mass in large objects, most 
collisions are catastrophic. Thus, the cratering algorithm has little impact on our results. 


Adopting the Bottke et al. (12010! ) treatment of m max ,d speeds up the collisional cascade. 
When the debris evolves more rapidly, swarms orbiting 10-30 M® planets cannot match the 
observed surface area of Fomalhaut b for ages of 100-400 Myr. However, more rapid evolution 
for swarms around 100-300 M® planets allows these systems to match the observations. 


5.2. Predictions for Fomalhaut b and Other Exoplanetary Systems 


Two aspects of our calculations allow tests from existing observations of Fomalhaut b 
or new observations of other debris disks. All collisional cascade models predict a mass 
loss rate from the production of part icles with sizes less than the size of the smallest stably 
orbiting particle (e.g.. iKobavashi fe Tanakall2010lh Most of these particle s should he close to 


the orbit of the planet around the central star (e.g., Kenyon et ah 2014). Numerical results 


for this mass loss rate at 100-400 Myr yield an expected surface brightness along the path 
of Fomalhaut b. For younger systems, the mass loss rate and the cross-sectional area of 


















































36 


the swarm are much larger. For sufficiently large Ad, satellite swarms are detectable around 
stars with ages of 1-10 Myr. 

To quantify our first prediction, we consider baseline models of satellite swarms orbiting 
10-100 M® planets. At 100-200 Myr, the mass loss rate in small particles is 0.6 — 3 x 
10 18 g yr -1 . For particle sizes of 100 /jrn, mass loss leaves behind a trail with a cross-sectional 
area of roughly 10 20 cm 2 every year. Fomalhaut b has an orbital period of roughly 1000 yr. 
Every orbit, mass loss produces a ring of material with a cross-sectional area comparable to 
the observed Ad of Fomalhaut b. 

The total surface area of this ring depends on the long-term evolution of small particles. 
If the particles have a velocity dispersion similar to their escape velocity from the planet, 
they have orbits with eccentricity e ~ 0.1 around Fomalhaut. Interactions with Fomalhaut 
b are probably rare. With inclinations i ~ e/2, the collision time for a single 100 /irn 
particle is roughly 1 Gyr for Ad & 10 23 and 1 Myr for Ad ~ 10 26 cm 2 . With roughly 1 Myr 
required to eject particles with Ad ~ 10 26 cm 2 , we envision an approximate steady-state 
where ejections of 100 /irn particles from Fomalhaut b roughly balance particles lost from 
destructive collisions. 


To estimate the surface brightness of this ring, we consider bound orbits along a ring 


with semimajor axis a ~ 120 AU, width 5a ~ 0.1 a, and e ~ 0.8 (e.g., iBeust et al 


Along this ring, there are roughly 10 4 resolution elements on HST images (Kenyon 


2014). 


et al. 


2014). With Ad ~ 10 26 cm 2 , each resolution element has a cross-sectional area of 10 22 cm 2 
in 100 /irn particles. 

If the orbit of Fomalhaut b is stable on Myr time scales, tracing a ring of dust alon g 


this orbit is challenging (e.g., ICurrie et al.l 120121; iGalicher et al.l 120131: iKalas et al.l 120131 1 . 


However, comparing the average surface brightness of coadded pixels along the orbit with 
similarly coadded pixels 20-30 AU away should yield a clear measure of the surface brightness 
along the ring and a strong test of the model. 

To make predictions for the brightness of satellite swarms orbiting any star, we consider 
f 0 the observed flux of the swarm relative to /* the observed flux from the central star 
(see also §2.3 of Keny on et aL- 2014). For a swarm with optical depth r and scattering 
efficiency Q s , f 0 /f* = QsAd/^na 2 . We set Ad = tA s = r7r(7]irHa p ) 2 with 771 = 0 . 2 . Defining 
Am = —2.5 log (/<>//*), the predicted contrast between the satellite swarm and the central 
star is 


Am « 15.83 - 2.5 log r - 5 log - 1.67 log 


M n 


10 Ma 


+ 1.67 log 


1M 0 


( 22 ) 


Although the contrast depends on r, 771 , M p , and M*, it is formally independent of the 
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semimajor axis of the planet and the distance to the star. However, swarms with r « 0.1-1 
have longer lifetimes at larger a; thus, observations are more likely to detect bright swarms 
at 100 AU than at 10 AU. 


Among nearby associations of young stars, the (3 Pic moving group provides the best 
testing ground fo r this prediction. With ~ 30 members having V « 4-9, d « 20 pc, and 
ages ~ 20 Myr ( Zuckerman et al. 2001 ; Mamaiek fe Belli 2014), this assoc iation is closer 
than the somewhat younger TW Hya (60 pc, 10 Mvr: [Schneider et al. 2012 ) and the some¬ 
what older Tuc-Hor (40 pc, 30 Myr; IZuckerman et al.l l201lh associat i ons. Roughly 30% 


of the members have luminous debris disks (e.g 


Rebull et al. 2008; Nilsson et al 


Riviere-Marichalar et al.ll2014T ) : [3 Pic contains at least one gas giant (ILagrange et al 


2009 


2010 ). 


From current planet detection statistics, > 50%-60% of stars with M+ < 1-2 have at 
least one planet with m p > 5-10 M 0 and a < 10-20 AU (e.g., Najita & Kenyon 2 0141 and 
references therein). If satellite swarms around massive planets are relatively common, there 
is a a reasonably high probability of finding at least one satellite swarm within the (3 Pic 
moving group. 

Although only one planetary mass companion has been det ected o rbiting members of 
the (3 Pic moving group, current detection limits are encouraging. 

Am ~ 9-10 mag at 075 in the broadband L filter; Biller et ah (2013) report Am 
at 1 - 2 " in the broadband H and narrow band CH 4 filters. As the sample sizes grow and 
the data acquisition/reduction techniques improve, it should be possible to constrain the 
frequency of luminous satellite swarms around the nearest young stars. 


Kasper et al.l (120071) derive 
14-15 


5.3. Predictions for Irregular Satellites in the Solar System 


Satellite swarm models for Fomalhaut b are based on the ense mble of irregular satellites 


orbiting the four gas giants in the solar system (e.g., Kennedy & Wyatt 2011). The ~ 160 


known satellites have radii r < 100-200 km and lie on ecc entric, high inclination orbits 


Brozovic et al. 

2011 ; 

Alexandersen et al. 

2012 ) 


Huong the largest objects with r fh 10 - 
100 km, the cumulative size distribution is shallow and reasonably close to a p ower law with 
n(> r) oc r~ Qc and q c ~ 1 (IJewitt fc Haghighipourl 120071: iBottke et al.l l2010l ) . For smaller 
objects, the size distribution may steepen. 


Using a set of coagulation calculations, Bottke et al. (120101 ) show that several Gyr of 
collisional evolution naturally produces satellite swarms with shallow size distributions. For 
model satellites with r fh 0.05-100 km orbiting Jupiter, the slope ranges from q c & 1.5 for 
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r < 5 km to q c < 1 at r = 5-100 km. Swarms with longer collision times orbiting Saturn 
and Uranus have somewhat steeper power laws. 

Our calculations confirm and extend these results. For 0.1-10 km satellites orbiting 
100-300 M® planets, size distributions at 1 Gyr are wavy power laws with q c ~ 1-2; 10- 
100 km objects have steeper power laws q c ~ 2-3. We extended several of these calculations 
to 5 Gyr; large objects then have q c ~ 1-2. 

Satellites with r < 0.1 km have wavy cumulative size distributions with a broad range of 
power law slopes f Figs. l9 MT6j) . Independent of various input parameters, very small particles 
with r < 1 cm have steep power laws with q c ps 4. Intermediate size particles with r ~ 1 cm 
to 100 m have flatter size distributions, q c ~ 0-2. These power laws are sensitive to the 
size of the smallest particles (Figs. fT2TfT3]i . Calculations with smaller particles have steeper 
power laws than calculations with larger particles. 


Although recent surveys detect several irregular satellites with r gg 1 km around Jupiter 


(e.g.j 

Brozovic et al. 2011; 

Alexander sen et al. 

2012 , 

Jacobson et al. 2012; 

Gomes-Junior et ah 

2015' 

, testing our predictions is challenging. 

With expected optical magnitudes > 28, ir- 


regular satellites with r < 0.1 km are too faint for any current and planned ground-based 
telescope. However, many of our calculations predict changes in the slope of the size distri¬ 
bution at 0.1-1 km. Extending the discovery space to this size range is feasible and would 
place interesting constraints on the coagulation models. 


6. CONCLUSIONS 

We describe results from a large suite of coagulation calculations for irregular satellite 
swarms orbiting 1-300 M® planets at a = 120 AU from a 1.9 M® central star. The calcula¬ 
tions follow the evolution of the size distribution for 10 /ini to 3000 km particles for 1 Gyr 
as a function of the initial mass of the swarm, the size of the smallest particle in the swarm, 
the initial size distribution of particles, the binding energy of the particles, and the method 
for distributing debris from a collision into smaller mass bins. 

Throughout the evolution, the largest satellites may grow or shrink. Growing satellites 
may scatter other satellites out of the planet’s Hill sphere or into tighter orbits around the 
planet. Among smaller satellites, the size distribution develops a characteristic shape with a 
steep power law at small sizes, a flat portion at intermediate sizes, and a shallow power law 
at larger sizes. The growth (shrinkage) of satellites and the time for the size distribution to 
develop a standard shape depend on the initial cloud mass, the initial size distribution, the 
initial r min and r max , and the binding energy of satellites. 
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In our baseline models, swarms orbiting 10-100 M® planets have cross-sectional areas 
at 100-400 Myr comparable to the observed cross-sectional area of Fomalhaut b. In these 
models, Xd = 0.01 and r min = 100 //m. Smaller Xd and larger r min allow swarms orbiting 
somewhat more massive planets to match observations of Fomalhaut b. Calculations with 
smaller r m j n require swarms around somewhat less massive planets. Changing the initial 
size distribution of satellites has little impact on these conclusions. Modifying the binding 
energy and the algorithm for distributing debris in smaller mass bins generally lowers the 
cross-sectional area, requiring swarms around more massive planets to match Fomalhaut b. 

Aside from discussing the impact of these calculations on our understanding of planet 
formation theory ( 1 15.II) . we derive predictions for (i) irregular satellites in the solar system 
and (ii) Fomalhaut b and satellites swarms in other exoplanetary systems. Identifying 0.1- 
1 km irregular satellites orbiting Jupiter would set interesting constraints on coagulation 
models. In Fomalhaut b, we predict a detectable trail of small particles within a few AU 
of the nominal orbit of the planet candidate. For exoplanetary systems with ages of 1- 
10 Myr, detectable satellite swarms orbiting 30-300 M® planets provide a way to estimate 
the frequency of sub-Jupiter mass planets at 50-150 AU around 1-2 M® stars. 


We acknowledge generous allotments of computer time on the NASA ‘discover’ cluster. 
Comments and suggestions from M. Geller, G. Kennedy, and an anonymous referee improved 
our discussion. Portions of this project were supported by the NASA Astrophysics Theory 
and Origins of Solar Systems programs through grant NNX10AF35G and the NASA Outer 
Planets Program through grant NNX11AM37G. 


A. Appendix 


To test the algorithms used in Orchestra , we compare numerical results with an alytic so- 


---— — y • — .. ~ 

lutions to the coagulation equation (Kenvon & Luu 

1998; Kenvon & Bromley 

2015) and pub- 

lished results from other investigators ( 

Kenvon & Bromlev 

2001 ; 

Brom 

ey & Kenvon 2006; 

Kenvon & Bromlev 2008; Bromlev & Kenvon 2011 

Kenvon & Bromlev 

2015 

). Here, we ex- 


amine how Orchestra performs for collisional cascades in spherical swarms of satellites orbit¬ 
ing a massive planet. 


The accuracy of all coagulation calculations de pends on the mass spacing parameter 


between adjacent mass b ins, 5k = m k+ i/m k fe.g.. IWetherilll Il990l: iKenvon A Luul Il998 


Kenyon A Bromley 2015). At the start of our calculations, we fix the typical mass m k 


and the boundaries 1/2 and m k+ i/ 2 of each mass bin. The initial average mass within 
each bin is fh k = M k /N k ; typically rh k ~ m k . As each calculation proceeds, collisions add 
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and remove mass from all bins; the average mass rrik and the average physical radius of 
particles fj. = (Sfhk/^'K p p ) 1 ^ then change with time. 

To illustrate how the evolution of satellite swarms changes with <5, we consider the 
baseline model described in the main text. An ensemble of satellites with Xd = 0.01, r max = 
100 km, and r m i n = 100 pm orbit a planet with M p = 10 M®. The initial size distribution 
is a power law with n(r ) oc r~ q and q = 3.5 between r m j n and r max . 

Fig. 1251 shows the time-variation of r max for five different values of 5. In the figure, all 
curves have the same general shape: a brief, ~ 10 5 — 10 6 yr period where r max is roughly 
constant, followed by a gradual decrease in r max with time. Tracks with larger 5 decline 
faster. 

Along each track, the decline consists of a gradual reduction in r max interspersed with 
occasional small jumps to larger r max and large jumps to smaller r max . In this example, 
collisions between equal mass objects with r ss r max increase the mass of the merged pair. 
These collisions produce jumps to larger r max . Cratering collisions - where somewhat smaller 
objects gradually chip away at the mass of larger objects - produce continuous mass loss 
from the largest objects. Thus, the average mass in the largest mass bin falls with time. 
Eventually, this mass falls below the mass boundary between adjacent bins (e.g., m^ < 
ixik— 1 / 2 ) • Objects in bin k are then placed into bin k — 1. Averaging the mass of the ‘old’ 
objects in bin k — 1 with the ‘new’ objects from bin k yields a new average mass fhk- 1 
which is smaller than the average mass of bin k. Thus, the size of the largest object jumps 
downward. Because the spacing of mass bins scales with 5, calculations with larger 5 have 
larger jumps than those with smaller 5. 

Although the mass loss rate from the grid is fairly insensitive to 5, the mass of the 
largest object clearly declines faster in calculations with larger 5. Cratering collisions are 
responsible for this difference. For all 5, these collisions are rare. Thus, only a few of the 
largest objects suffer substantial mass loss from cratering collisions every time step. When 
<5 is small (1.05-1.10), these objects are placed into the next smallest mass bin; the average 
mass of the remaining objects in the mass bin is unchanged. When 5 is large (1.41-2.00), 
the amount of mass loss is not sufficient to place objects into the next smallest mass bin; 
the average mass of all objects in the bin then decreases. As a result, the average mass of 
the largest objects declines faster when 5 = 2 than when 5 = 1.05. 

Despite this difference, other aspects of the evolution are fairly insensitive to 5. Fig. [221 
shows snapshots of the relative size distributions at 100 Myr. Each curve follows a standard 
pattern, with a steep power law at 0.1 mm to 10 cm, a minimum at ~ 30 cm, a rise from 1 m 
to 100 m, a wavy pattern from 100 m to 50 km, and then an abrupt decline at the largest 
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sizes. When <5 = 2, the fluctuations about a reference model with 5 = 1.05 are large. For 
other 5, deviations from the reference model are minimal. 

Fig. EDI illustrates the evolution of the relative surface area for two different baseline 
models as a function of 5. When M p = 10 M e , the relative area declines from roughly 10 2 
at 100-1000 yr to roughly 0.1 at 1 Gyr. When 5 = 1.05, the decline is smooth, with a minor 
change in slope at roughly 10 4 yr. Adjustments from the initial power law to the non-power 
law equilibrium size distribution (e.g., Fig. [2D]) cause this change in slope. Evolution with 5 
= 2 is more ragged, with modest fluctuations relative to the reference model with 5 = 1.05. 
As 5 declines, the evolution of the relative area follows the reference model more closely. 

When M p = 100 M e , the initial relative area for models with Xd = 0.01 is a factor 
of ten larger. The long term evolution is similar: a slow decline with an inflection point 
around 10 4 yr. Once again, the evolution of the relative area is somewhat more ragged in 
calculations with 5 = 2 than in calculations with smaller 5. 

Although there are clear differences in the evolution as a function of 5, the ability of 
an initial set of model parameters to match the observations rarely depends on 5. For these 
two examples, all of the 100 M p calculations pass through the target box for Fomalhaut b. 
All of the 10 M p calculations graze the lower edge of the target box. 

For swarms of satellites in a spherical shell around a massive planet, calculations with 
5 < 1.2 yield a better understanding of the long term evolution of r max and the size distribu¬ 
tion. Evolution of the total mass and relative surface area are fairly insensitive to 5. Single 
annulus calculations with 5 = 1.2 run quickly, with execution times of 25 cpu hours for 1 Gyr 
evolution times using 10 /ini to 1000 km particles. Thus, we perform most calculations with 
5 = 1.2 and use occasional calculations with smaller 5 to verify interesting features of the 
evolution. 
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Fig. 1.— Time evolution of the relative surface area for the analytic model with a = 1 
described in the text. The legend in the lower left indicates M p the mass of the planet in M® 
for each solid curve. The legend in the upper right indicates Xd the mass of circumplanetary 
material relative to the mass of the planet and r max the radius (in km) of the largest object 
in the swarm. Eq. [8] sets values for r m ; n the radius (in //in) of the smallest particle in the 
swarm. The grey shaded box indicates the locus of allowed points for Fomalhaut b, using 
the surface area derived from the measured brightness, the age of Fomalhaut, and la errors. 
For the adopted combination of Xd and r maX) clouds orbiting planets with M p ~ 30-100 M® 
match the observations. Models with M p = 10 M® and 300 M® barely miss the shaded box; 
those with M p = 1 M® have too little surface area at all times. 
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Fig. 2.— Time evolution of the radius of the largest object derived from coagulation calcu¬ 
lations of circumplanetary clouds of particles with the nominal fragmentation parameters, 
Xd = 0.01, r min = 100 fjL m, m/ j0 = 0.2, and bi = 0.0. The legend in the upper left indicates 
M p for each solid curve. For calculations with massive planets, the size of the largest object 
smoothly declines with time. When the mass of the planet is smaller, the largest objects 
sweep up small particles and grow into much larger objects. In calculations where r max is 
400 km (50 km), the largest objects are more (less) likely to grow with time. 
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Fig. 3.— As in Fig. [2] for calculations starting with r max = 50 km and 200 km. When M p 
is large, the collisional cascade gradually destroys the largest objects. Around lower mass 
planets, the largest objects grow with time. 
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Fig. 4.— As in Figs. [2H3] for calculations with M p = 300 M® and x<i as indicated in the 
legend. When the mass of the cloud is smaller, the radius of the largest object changes more 
slowly. 
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Fig. 5.— As in Figs. EH3] for calculations with M p = 10 M e and r min = 10 /mi (orange 
curves), r min = 100 /mi (green curves), and r min = 1000 /mi (1 mm; blue curves). The 
growth of large objects is fairly independent of the size of the smallest particles in the grid. 
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Fig. 6.— As in Figs. EH3]for calculations with M p = 10 M®, 0 — 0.2, and either bi = 0.0 

(blue curves) or 5; = 1.0 (orange curves). The growth of large objects is fairly independent of 
the exponent in the relation between the mass of the largest object and the collision energy. 
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Fig. 7.— As in Figs. [2H3] for calculations with M p = 10 M e and different Q* D . The legend 
indicates the value of Q* D relative to the nominal fragmentation parameters listed in the 
main text. 
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Fig. 8.— As in Figs. [2H3] for calculations with M p = 10 M e and different initial size 
distributions. The legend indicates whether the calculation starts with a mono-disperse 
set of satellites (no sd) or a power law size distribution. When the initial population is 
mono-disperse, satellites grow faster. 
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Fig. 9.— Snapshots evolution of the relative cumulative size distribution for calculations 
with M p = 10 M 0 , Xd = 0.01, r max = 100 km, r min = 100 pm, m; )0 = 0.2, and b t = 0.0. 
The legend in the upper right indicates the evolution time in Myr. Within roughly 1 Myr, 
collisions produce several distinct features in the relative size distribution: (i) a steep rise at 
the largest sizes (r > 50 km), (ii) a shallower rise which approximately follows the original 
power law, (iii) a sharp drop at intermediate sizes (r fa 30 cm to 0.1 km), and (iv) a steep 
rise at the smallest sizes (r fa 0.1 mm to 30 cm). 
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Fig. 10.— As in Fig. [9]for a range of masses for the central planet at t — 100 Myr. The legend 
in the upper right indicates the mass of the planet in M 0 . For all M P) the size distribution is 
very steep for particle sizes r « 0.1 mm to 1-30 cm. The depth of the minimum at 10-30 cm 
grows with the mass of the planet. For r > 0.1 km, fluctuations about the original power 
law grow with the mass of the planet. 
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Fig. 11.— As in Fig. ITOlfor M p = 10 M 0 and various initial r max as indicated in the legend. 
Aside from differences at r ^ r max , the relative size distribution at 100 Myr is independent 
of initial r max . 






Fig. 12.— As in Fig. |TT] for r max = 100 km and various r min as indicated in the legend. 
When r min is smaller, the relative size distribution is closer to a single power law for r > 
1 cm and has a smaller deficit of particles at 10-1000 cm. 
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Fig. 13.— As in Fig. [TO] for M p = 10 M® and various combinations of r maX) m.L,o, and b L as 
indicated in the legend. For r > 10 3 cm, the relative size distribution is fairly independent 
of &£. Among smaller particles, calculations with bj J — 1 yield a shallower size distribution 
than those with b^ = 0. 
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Fig. 14. — As in Fig. [TT| for r max = 100 km, r min = 100 /iin, and different values for Q* n . The 
legend indicates the value of Q* D relative to the nominal fragmentation parameters listed in 
the main text. 






61 



Radius (cm) 


Fig. 15.— As in Fig. [9] for calculations starting with a mono-disperse set of particles. The 
legend indicates the evolution time in Myr. At early times, collisions among 100 km objects 
produce a small amount of debris populating the small size end of the size distribution. After 
10 Myr, debris from these collisions and the debris from collisions of smaller objects yields 
a smooth equilibrium size distribution from 10 /im to roughly 100 km. 
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Fig. 16.— Comparison of model relative size distributions for calculations starting from an 
initial power law size distribution of particles (‘sd’) and a mono-disperse set of particles (‘no 
sd’) at 10 Myr (upper set of curves) and at 1 Gyr (lower set of curves). Aside from minor 
deviations at the largest sizes, the two sets of calculations yield identical size distributions. 
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Fig. 17.— Time evolution of the relative surface area derived from coagulation calculations 
of circumplanetary clouds of particles with the nominal fragmentation parameters, o = 
0.2, and b L = 0.0. The legend in the upper right indicates the initial xj, r max (in km), and 
r mm (in /iin) . The legend in the lower left indicates the mass of the central planet. When 
the largest object does not grow (M p = 100, 300 M®), the relative surface area declines 
smoothly with time. At late times in simulations with growing satellites, the surface area 
fluctuates about a gradual decline. 
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Fig. 18.— As in Fig. [17] for the M p , Xd, and r min indicated in the upper right corner for 
the range of r max (in km) indicated in the lower left corner. The slow decline of the relative 
surface area is smoother for smaller r max . 
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Fig. 19.— As in Fig. [17] for the M p , r max , and r min indicated in the upper right corner 
for various Xd as indicated in the lower left corner. Lower mass clouds have smaller relative 
surface areas. 
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Fig. 20.— As in Fig. [T9] for M p = 100 M®. 
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Fig. 21.— As in Fig. [17] for the M pi and r max summarized in the upper right corner 
for the r m i n listed in the lower left corner. Calculations with smaller r m i n have more small 
particles and larger relative surface area. 
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Fig. 22.— As in Fig. [21] for the Xd, r max , and r min listed in the upper right corner and 
various M p and as listed in the lower right corner. Calculations with = 1 have smaller 
relative surface area than those with = 0. 
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Fig. 23.— As in Fig. [17] for the M p , xd, r max , and r min listed in the upper right corner and 
various Q* D . The legend in the lower left corner indicates the value of Q* D relative to the 
nominal value. 
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Fig. 24.— As in Fig. [17] for calculations with (sd) and without (no sd) an initial power 
law size distribution of solid particles. The legend indicates M p , r max , and the initial size 
distribution for each model curve. 
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Fig. 25.— Evolution of semimajor axis for massive (n-body) satellites orbiting a 1 M® 
planet. The dashed lines indicate the extent of the satellite swarm within the coagulation 
code. Tracks for individual n-bodies are coded by color. At 70-100 Myr, the coagulation 
code promotes three objects into the n-body code. After several minor encounters, strong 
interactions lead to a single object on an e = 0.6 orbit at smaller a and the ejection of two 
n-bodies. Somewhat later (~ 200 Myr), interactions between a second pair of n-bodies leads 
to a second set of ejections and a modest contraction of the orbit of the original n-body. At 
late times, promotion of a sixth n-body leaves the system with a single relatively low mass 
object on a fairly circular orbit at 0.20 AU and a more massive satellite on an eccentric orbit 
at 0.05 AU. 
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Fig. 26.— Evolution of the mass of satellite swarms relative to their initial mass for the 
baseline numerical model (’num’) and analytic models (’an’) with three values of a as listed 
in the legend. As indicated in the upper right corner, the baseline model has M p = 10 M®, 
Xd = 0.01, r max = 100 km, and r min = 100 /irn. The evolution of an analytic model with 
a ~ 1 provides a reasonable match to the numerical model. 
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Fig. 27.— Grid of outcomes for satellites swarms in Fomalhaut b. For initial relative cloud 
mass Xd = 0.001 and 0.01 (as indicated above the first column of points), green (red) points 
indicate models which match (do not match) the observed surface area of Fomalhaut b at 
100-400 Myr. 
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Fig. 28.— Time evolution of the radius of the largest object for different mass spacing 
factors 5 as listed in the legend. The satellite swarm - x,i — 0.01, r max = 100 km, r min = 
100 /mi - orbits a planet with M p = 10 M®. 
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Fig. 29.— Snapshot of the size distribution for the baseline models in Fig. [28] at 100 Myr. 
Solutions for 6 = 2 have more waviness than those with smaller 5. Solutions for 6 = 1.05-1.4 
are nearly identical. 
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Fig. 30.— Time evolution of the cross-sectional area for two baseline models as a function 
of S (as listed in the legend). For the upper (lower) set of curves, M p = 100 M® (10 M®); 
in both, Xd = 0.01, r max = 100 km, and r mm = 100 /mi. The surface area is not a strong 
function of 5. 








